I am working with a dataset that includes raster data that I'm combining with maps from rnaturalearth. Most of what I'm doing is manipulating the raster data using terra, to do things like aggregation, reprojection, and then various raster operations using app (for things like taking the mean of multiple layers, etc).
For various visualization purposes, I sometimes want to crop my rasters and display regional data using ggplot. Cropping the raster itself is easy using crop(), but cropping the sf object is a little harder. The solution I found in older tutorials now throws an error:
library(sf)
library(rnaturalearth)
world <- ne_countries(scale = "small", returnclass = "sf")
europe <- st_crop(world, xmin = -20, xmax = 45,
ymin = 30, ymax = 73)
Error in wk_handle.wk_wkb(wkb, s2_geography_writer(oriented = oriented, :
Loop 0 is not valid: Edge 0 crosses edge 78
Apparently this is caused by the switch from flat space to spherical geometry, and can be solved by turning off sf_use_s2():
sf_use_s2(FALSE)
europe <- st_crop(world, xmin = -20, xmax = 45,
ymin = 30, ymax = 73)
Spherical geometry (s2) switched off
although coordinates are longitude/latitude, st_intersection assumes that they are planar
Warning: attribute variables are assumed to be spatially constant throughout all geometries
This vignette goes into a lot of detail about spherical geometry in sf, but I'm not sure if I'm catching all the nuances, and I'm afraid of inadvertantly doing something to mess up the geometry without realizing it. In general, I'm working with lat/lon (WGS84) data, and I do all my analysis with lat/lon, then project to mollweide using coord_sf() for visualization purposes only. I'm largely working with global data, which I occassionaly visualize at more local scales. So my question is:
Is there any problem in turning off sf_use_s2() for the purposes of cropping sf objects? Should I turn it back on after performing this operation?
This warning
although coordinates are longitude/latitude, st_intersection assumes that they are planar
Can normally be ignored when using the data for display (making a map). And more generally it may not be very important, as long as the distance between nodes (vertices) in your lines/polygons is not very large. You can increase the minimum node distance with terra::densify (perhaps there is an equivalent in sf).
The message could be read as suggesting that the problem would go away if you transform the data to a planar projection, but that is not the case. The message would go away, but the possible imprecision would not be fixed.