Find maximum extent out of list of shapefiles (in projected coordinate system) in R

Question

I have a script where I crop a list of rasters by the extent of one of many shapefiles. So far all of by shapefiles have been such that I could manually select the one with the largest extent (i.e. all others fit within the largest one). But I'm trying to make my script automated (and robust) so that it generates the largest extent considering a list of shapefiles that may overlap (i.e. one has greatest x_max, another greatest x_min).

Here's a (mostly) reproducible example (even though birds is within cities, it's a start). The output of my code would ideally generate something with the same as the extent of cities (a better example would be one with mixed values, but this is the best I could find).

library(rgdal)
library(raster)

dsn <- system.file('vectors', package = 'rgdal')[1]
cities <- readOGR(dsn = dsn, layer = 'cities')
birds <- readOGR(dsn=dsn, layer = "trin_inca_pl03")

shapefile_vector <- c(cities, birds)
shapefile_plots <- lapply(X = shapefile_vector, FUN = plot)

extent(cities)
extent(birds)

find_extent <- function(shapefile){
  as.matrix(extent(shapefile))
}

shapefile_extents <- lapply(X = shapefile_vector, FUN = find_extent)
shapefile_extents 

I can't figure out how to go from a list of matrices (shapefile_extents) to a list of absolute max coordinates to use in something like (taken from Crop a raster file in R):

e <- as(extent(-16, -7.25, 4, 12.75), 'SpatialPolygons')
crs(e) <- "+proj=longlat +datum=WGS84 +no_defs"
r <- crop(worldpopcount, e)

Answer 1

Given a list of shapey objects:

> shapes = list(cities, birds)

Compute the list of extents of those shapey objects

> shape_extents = lapply(shapes, raster::extent)

and then call raster::merge on them to get the merged extent:

> do.call(raster::merge, shape_extents)
class      : Extent 
xmin       : -165.27 
xmax       : 177.1302 
ymin       : -53.15 
ymax       : 78.2 

I think that saves you having to do the matrix conversion and then manipulate that.

Note this is flexible, in that it works on any object with an extent method (which includes rasters, sp and sf objects, be they points, lines, polygons or anything else), and minimal, in that it doesn't do any geometry reorganisation or conversion that might slow things down.

Answer 2

The raster::extent() can take SpatialObject as the input (not only (xmin, xmax, ymin, ymax)).

If the input geometries are the same type, extent(union(cities, birds)) should work. However, as the cities are points and the birds are polygons in the provided example, we need a workaround (I would suggest rgeos::gUnion()).

# The CRS have to be the same, but unfortunately birds' CRS is NA.
crs(cities)  # WGS84
crs(birds)   # NA
crs(birds) <- "+proj=longlat +datum=WGS84 +no_defs"

library(rgeos)
e <- as(extent(rgeos::gUnion(cities, birds)), 'SpatialPolygons')
extent(e)
# class      : Extent 
# xmin       : -165.27 
# xmax       : 177.1302 
# ymin       : -53.15 
# ymax       : 78.2 

Next step is the same as the one you have linked to.

crs(e) <- "+proj=longlat +datum=WGS84 +no_defs"
r <- crop(data_to_be_cropped, e)

Answer 3

So if you have some extent objects from your shapefiles (ext1,ext2,ext3...) then you get full extent:

spoly <-  do.call(bind, sapply(c(ext1,ext2,ext3), FUN = function(x){as(x, 'SpatialPolygons')})) 

full_extent <- extent(bbox(spoly))