I have multipolygons that span across the globe. I was wondering if it was possible to re-project each row, i.e. each country/ ISO3 to a local projection, using a central grouping mechanism.
For example, my polygons within the a particular region say SE Asia is then re-projected to a local projection for that area. Preferably to a LAEA type projection as I want to create buffers around these polygons for further analyses.
Here is the structure of my data and the projection I have used for both multipolygons and multipoints.
Simple feature collection with 6 features and 2 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: -11959630 ymin: -4449011 xmax: 5370704 ymax: 10196680
epsg (SRID): NA
proj4string: +proj=laea +lat_0=45.5 +lon_0=-114.125 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs
PARENT_ISO ISO3 geom
1 ABNJ ABNJ MULTIPOLYGON (((-7616720 -4...
2 ARE ARE MULTIPOLYGON (((2426453 101...
3 ATG ATG MULTIPOLYGON (((5364022 -14...
4 AUS AUS MULTIPOLYGON (((-10014124 -...
5 AUS CCK MULTIPOLYGON (((-8801631 79...
6 AUS CXR MULTIPOLYGON (((-9722602 63...
> head(c_pt)
Simple feature collection with 4 features and 2 fields
geometry type: MULTIPOINT
dimension: XY
bbox: xmin: -409654.4 ymin: -3100315 xmax: 3235020 ymax: -1443240
epsg (SRID): NA
proj4string: +proj=laea +lat_0=45.5 +lon_0=-114.125 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs
PARENT_ISO ISO3 geom
1 BLZ BLZ MULTIPOINT (2757788 -286421...
2 GTM GTM MULTIPOINT (2773602 -287236...
3 HND HND MULTIPOINT (2883494 -283266...
4 MEX MEX MULTIPOINT (-409654.4 -1831...
Creating buffers from this projection seemed to create a lot of distortion, as can be seen by Madagascar.
I do realise I could split the data by country and reproject each new country file, but was wondering if there was a way that was possible whilst keeping all the data together?
Or is this just not the way to handle this kind of situation?
Example of my workflow
#Reproject coral reefs
c_pt <- st_transform(c_pt, crs = "+proj=laea +lat_0=45.5 +lon_0=-114.125 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs")
c_py <- st_transform(c_py, crs = "+proj=laea +lat_0=45.5 +lon_0=-114.125 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs")
#Create buffer of 100km around coral reefs
buffer_100_py <- st_buffer(c_py, dist = 100000)
buffer_100_pt <- st_buffer(c_pt, dist = 100000)
#this is gives overlapping portion of polygons
overlap <- st_intersection(buffer_100_py, buffer_100_pt)
#this is to delete the overlapping polygons
diffPoly <- st_difference(buffer_100_pt, st_union(overlap))
#joining the cleaned point buffer data with the polygon data
buffer_100_all<-rbind(diffPoly, buffer_100_py)
#Now to clean up the polygons to remove overlap by union by feature - there are still overlaps of polygons,
#But this should be fine for the population extraction
buffer_100_all_one<-ms_dissolve(buffer_100_all, field = "ISO3")
UPDATE:
Using the row-wise projection example provided by @mdsumner, I was able to create 100 km buffers with much less/no visible distortion.
My workflow was as follows
#Functions provided by @mnsumner
local_proj <- function(lonlat) {
sprintf("+proj=laea +lon_0=%f +lat_0=%f +datum=WGS84", lonlat[1], lonlat[2])
}
local_reproj <- function(x) {
cc <- sf::st_coordinates(st_centroid(x))
sf::st_transform(x, local_proj(cc))
}
## easiest way to row-wise is split; here applying local projections to each feature
c_py_row <- lapply(split(c_py, 1:nrow(c_py)), local_reproj)
##create buffer function to apply to each row-wise
##Function for 100km here
buffer_100 <- function(x) {st_buffer(x, dist = 100000)
}
#buffering each row by 100km
buffer_100_py_row <- lapply(c_p_row, buffer_100)
#function to re-project each row to EPSG 4326
world_crs <- function(x) {st_transform(x, crs = 4326)
}
##re-projecting each row to same CRS
buffer_100_py_row_2<-lapply(buffer_100_py_row, world_crs)
#binding rows back into single dataframe
buffer_100_py_row_3<-do.call(rbind, buffer_100_py_row_2)
Similar to the country wise plot in the GitHub link, here are my polygons from the row-wise split.
