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_bbox
and st_transform
). In contrast, my
current experience of trying to use terra
and 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 raster
or terra
. Some colleagues
told me they’ve had unresolvable issues with stars
and reverted to
raster
. I don’t know the details though.
There is a ‘babel’ for stars to raster function names though.
Data input and plots are quite similar to sf
:
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 sf
:
#create an extent object
ext2 <- st_bbox(r)
#constrain it in x direction
ext2[1] <- 120
ext2[3] <- 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
So terra
was about 8.5 times faster than raster
and stars
was
about 3 times faster than terra
(so about 26 times faster than
raster
).
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
find out).
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.
Designed by Chris Brown. Source on Github