How to crop non-projected EPGS:4326 polygons along straight lines?

Viewed 102

I am trying to allocate spatial polygons to unevenly-sized boxes defined by decimal degree coordinates with the sf package:

library(sf)
#> Linking to GEOS 3.8.1, GDAL 3.2.1, PROJ 7.2.1
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(spData)
#> To access larger datasets in this package, install the spDataLarge
#> package with: `install.packages('spDataLarge',
#> repos='https://nowosad.github.io/drat/', type='source')`

x <- world 

plot(st_geometry(x), xlim = c(60,180), ylim = c(0, 60))

matrix(c(60, 0, 180, 0, 180, 30, 60, 30, 60, 0), byrow = TRUE, ncol = 2) %>%
  list() %>% st_polygon() %>% st_sfc() %>% st_set_crs(4326) %>% 
  plot(add = T, border = "red", col = NA, lwd = 3)

matrix(c(60, 30, 120, 30, 120, 60, 60, 60, 60, 30), byrow = TRUE, ncol = 2) %>%
  list() %>% st_polygon() %>% st_sfc() %>% st_set_crs(4326) %>% 
  plot(add = T, border = "blue", col = NA, lwd = 3)

I want to crop/split the contents of the blue and red boxes into separate entities. When I use the st_crop function, I notice that the cropping is done along curved lines coinciding with the edges of the underlying polygons of x:

crop1 <- st_crop(x, xmin = 60, xmax = 180, ymin = 0, ymax = 30)
#> Warning: attribute variables are assumed to be spatially constant throughout all
#> geometries
crop2 <- st_crop(x, xmin = 60, xmax = 120, ymin = 30, ymax = 60)
#> Warning: attribute variables are assumed to be spatially constant throughout all
#> geometries

plot(st_geometry(x), xlim = c(60,180), ylim = c(0, 60))
plot(st_geometry(crop1), border = "red", col = NA, add = TRUE)
plot(st_geometry(crop2), border = "blue", col = NA, add = TRUE)

Created on 2021-09-02 by the reprex package (v2.0.1)

It seems that the coordinates are not handled planar as stated here, but perhaps projected before cropping (?). Also, one could reduce the curvature by adding nodes to the underlying polygons before cropping. While this may be required when further projecting the cropped object, it seems unnecessarily complicated way as opposed to just doing the cropping by assuming lon/lat values as cartesian coordinates.

Is there (an easy) way to crop the decimal degree polygons assuming cartesian coordinates (i.e. along straight lines) in sf? I am mostly concerned by the resulting overlap of polygons when they are uneven in size.

EDIT

sp and rgeos do not seem to have this problem:

``` r
library(sf)
#> Linking to GEOS 3.8.1, GDAL 3.2.1, PROJ 7.2.1
library(sp)
library(rgeos)
#> rgeos version: 0.5-7, (SVN revision (unknown))
#>  GEOS runtime version: 3.9.1-CAPI-1.14.2 
#>  GEOS using OverlayNG
#>  Linking to sp version: 1.4-5 
#>  Polygon checking: TRUE
library(spData)
#> To access larger datasets in this package, install the spDataLarge
#> package with: `install.packages('spDataLarge',
#> repos='https://nowosad.github.io/drat/', type='source')`

y <- as_Spatial(world)
rect1 <- matrix(c(60, 0, 180, 0, 180, 30, 60, 30, 60, 0), byrow = TRUE, ncol = 2) %>%
  list() %>% st_polygon() %>% st_sfc() %>% st_set_crs(4326) %>% as_Spatial()
rect2 <- matrix(c(60, 30, 120, 30, 120, 60, 60, 60, 60, 30), byrow = TRUE, ncol = 2) %>%
  list() %>% st_polygon() %>% st_sfc() %>% st_set_crs(4326) %>% as_Spatial()

sp::plot(y, xlim = c(60,180), ylim = c(0, 60))
sp::plot(rgeos::gIntersection(y, rect1, byid = TRUE), add = TRUE, border = "red", col = NA)
sp::plot(rgeos::gIntersection(y, rect2, byid = TRUE), add = TRUE, border = "blue", col = NA)

Created on 2021-09-02 by the reprex package (v2.0.1)

0 Answers
Related