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
     )
spatial-subsets_9_0.png 
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")
spatial-subsets_13_0.png