The other day I wrote a comparison of raster and terra packages. The other newish raster package on the block is stars. So I go through the same code here.
If you are starting out with stars, read their manual first.
But first some key differences.
stars is set-up to work with spatio-temporal arrays (terra also has better capacity to work with multilayer rasters than raster did). stars deals with more complex data types than either raster or terra, for instance, stars can handle rotated grids.
stars works will with
sf, and many
sf functions have methods for
stars objects (like
st_transform). In contrast, my
current experience of trying to use
sf is a bit more
awkward. I had to do lots of class conversions of
terra objects to
sf objects and back again.
stars isn’t as well documented as
terra. Some colleagues
told me they’ve had unresolvable issues with
stars and reverted to
raster. I don’t know the details though.
Data input and plots are quite similar to
library(stars) r <- read_stars("data-for-course/spatial-data/MeanAVHRRSST.grd") plot(r)
st_bbox(r) ## xmin ymin xmax ymax ## 82.50 -72.25 181.25 -9.75
Now let’s crop and reproject it, it looks almost identical to
reprojecting a polygon in
#create an extent object ext2 <- st_bbox(r) #constrain it in x direction ext2 <- 120 ext2 <- 170 r2 <- st_crop(r, ext2) r3 <- st_transform(r2, "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs +towgs84=0,0,0") plot(r3)
Let’s compare the three packages for speed at projecting data, an often time-consuming task.
library(microbenchmark) r_terra <- terra::rast("data-for-course/spatial-data/MeanAVHRRSST.grd") r_raster <- raster::raster("data-for-course/spatial-data/MeanAVHRRSST.grd") robin_proj <- "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs +towgs84=0,0,0" tout <- microbenchmark( raster::projectRaster(r_raster, crs=robin_proj), terra::project(r_terra, robin_proj), st_transform(r, crs = robin_proj), times = 10 ) tout ## Unit: milliseconds ## expr min lq ## raster::projectRaster(r_raster, crs = robin_proj) 330.49234 348.60983 ## terra::project(r_terra, robin_proj) 66.47145 67.78466 ## st_transform(r, crs = robin_proj) 19.03904 19.91101 ## mean median uq max neval ## 506.93616 505.29937 521.27790 874.32388 10 ## 69.90088 68.87271 69.09739 80.66674 10 ## 23.38595 20.64656 22.88241 42.57163 10
terra was about 8.5 times faster than
about 3 times faster than
terra (so about 26 times faster than
The answer here obviously depends on what packages you want to use. A key one for me is tmap for mapping:
library(tmap) tm_shape(r) + tm_raster()
Apparently you can
coerce stars objects into raster objects, but I
couldn’t figure out easily how to do that (with the 10 minutes I had to
Its a confusing world of raster packages now, with raster, terra and stars all viable options. A decision about your ‘go to’ package needs more exploration than my brief blog and will depend on what types of operations you are doing.
The real test will come when analyzing large and complex datasets. The comparison I’ve made here is very shallow, it doesn’t explore the full functionality of terra or stars, which is much deeper than raster. So I’m interested to hear from folk how they go with more sophisticated applications.