I have a following use case to cut/diff the polygons that spread outside the boundary. I am trying to migrate from rgdal package to sf packages as they are expiring. The code below works fine.
WorldCountry <- sf::st_read("data/modified.countries.geo.json")
countries <- c("Canada", "Mexico")
data_Map <- WorldCountry[WorldCountry$name %in% countries, ]
data_Map <- as_Spatial(data_Map)
for (i in seq(length(spa_polys), 3)) {
spa_polys@proj4string <- data_Map@proj4string
geo_diff <- gDifference(spa_polys[i], data_Map)
iso@polygons[[i]] <- geo_diff@polygons[[1]]
iso@polygons[[i]]@ID = as.character(i)
rm(geo_diff)
}
It generates the image with cut polygons. However, while changing the code to use SF packages, it does not correctly calculates the difference for both st_intersection and st_difference.
WorldCountry <- sf::st_read("data/modified.countries.geo.json")
countries <- c("Canada", "Mexico")
data_Map <- WorldCountry[WorldCountry$name %in% countries, ]
data_Map <- as_Spatial(data_Map)
data_Map_sf <- st_as_sf(data_Map, coords = c("longitude", "latitude"), crs = 4326)
for (i in seq(length(spa_polys), 1)) {
spa_poly_sf <- st_as_sf(spa_polys[i])
st_crs(spa_poly_sf) <- st_crs(data_Map_sf)
if (!st_is_longlat(spa_poly_sf)) {
spa_poly_sf <- st_transform(spa_poly_sf, st_crs(data_Map_sf))
}
geo_cut <- st_intersection(spa_poly_sf, data_Map_sf)
iso$geometry[[i]] <- geo_cut
iso$ID[i] <- as.character(i)
rm(geo_cut)
}
As shown in the figure, the polygons are not cut.
Any suggestions would be greatly helpful.
I have tried both methods using the intersection and difference. I am expecting it to work like it worked with rgdal.