Subsetting spatial feature objects¶
The following makes use of the sf package (in a version prior to 1.0). You may need to install it from
CRAN using the code
install.packages("http://cran.r-project.org/src/contrib/Archive/sf/sf_0.9-7.tar.gz",repos=NULL)
if you want to run this on your computer.
library(sf)
Linking to GEOS 3.9.0, GDAL 3.2.2, PROJ 7.2.1
The files “south-america-1990.RData” and “ged191.RData” were created in previous examples.
load("south-america-1990.RData")
load("ged191.RData")
# This fails due to different coordinate reference systems
Colombia.conflicts <- ged191[Colombia,]
Error in st_geos_binop("intersects", x, y, sparse = sparse, prepared = prepared, :
st_crs(x) == st_crs(y) is not TRUE
st_crs(Colombia)
Coordinate Reference System:
User input: +proj=longlat +ellps=WGS84
wkt:
GEOGCRS["unknown",
DATUM["Unknown based on WGS84 ellipsoid",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1],
ID["EPSG",7030]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8901]],
CS[ellipsoidal,2],
AXIS["longitude",east,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433,
ID["EPSG",9122]]],
AXIS["latitude",north,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433,
ID["EPSG",9122]]]]
st_crs(ged191)
Coordinate Reference System:
User input: EPSG:4326
wkt:
GEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
USAGE[
SCOPE["Horizontal component of 3D system."],
AREA["World."],
BBOX[-90,-180,90,180]],
ID["EPSG",4326]]
# Here we make sure that the coordinate reference systems match
ged191_ellips <- st_transform(ged191,st_crs(Colombia))
Colombia.conflicts <- ged191_ellips[Colombia,]
although coordinates are longitude/latitude, st_intersects assumes that they are planar
although coordinates are longitude/latitude, st_intersects assumes that they are planar
mypal <- function(n)gray.colors(n,start=.2,end=.9,alpha=.5)
plot(st_geometry(Colombia))
plot(Colombia.conflicts["deaths_civilians"],
add=TRUE,pch=19,cex=.2,
pal=mypal,
nbreaks=30
)

st_circ <- function(x,dist.km){
dist.degr <- 360*dist.km/40007.863
st_buffer(st_geometry(x),dist=dist.degr)
}
Bogota.200km <- st_circ(Bogota,dist.km=200)
Warning in st_buffer.sfc(st_geometry(x), dist = dist.degr):
st_buffer does not correctly buffer longitude/latitude data
dist is assumed to be in decimal degrees (arc_degrees).
Bogota.conflicts <- ged191_ellips[Bogota.200km,]
although coordinates are longitude/latitude, st_intersects assumes that they are planar
although coordinates are longitude/latitude, st_intersects assumes that they are planar
# The following seems to work only for sf version < 1.0 :(
plot(st_geometry(Colombia))
plot(st_geometry(Colombia.conflicts),
add=TRUE,pch=1,cex=.3,col="gray80"
)
plot(Bogota.200km,lty=3,add=TRUE)
plot(st_geometry(Bogota.conflicts),
add=TRUE,pch=19,cex=.3,col="black")

- R file: spatial-subsets.R
- Rmarkdown file: spatial-subsets.Rmd
- Jupyter notebook file: spatial-subsets.ipynb
- Interactive version of the Jupyter notebook (shuts down after 60s):
- Interactive version of the Jupyter notebook (sign in required):