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.

enter image description here

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")


Using the row-wise projection example provided by @mdsumner, I was able to create 100 km buffers with much less/no visible distortion.

enter image description here

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.

enter image description here

  • This is all possible, but there's at least a couple of problems here - first, how to choose a local projection given lat-long, and second how to store data reprojected to different CRS. Once the first problem is solved you can add the CRS as a column to the data, and then when you want to do stuff with that CRS work row-wise, apply the transformation, do the thing, and then transform back if you want. Its all possible. You need to code it though.
    – Spacedman
    Sep 11, 2019 at 13:25
  • A short example: gist.github.com/mdsumner/5af323f455d839b80c41bb043f5b2068 I thought this a perfectly good question and shouldn't be closed. There's no help out there for this, and why not contribute some.
    – mdsumner
    Sep 17, 2019 at 18:55
  • @Spacedman thank you for confirming the possibilities.
    – Asw
    Sep 23, 2019 at 13:22
  • @mdsumner, thank you so much for your example! I am currently running this on my dataset. I did wonder if you had any advice on combining the output of "listof" back into a data frame as they have differing CRS.
    – Asw
    Sep 23, 2019 at 13:24
  • well, you can do that with a data frame and a list column, but it's probably not a great idea - there's a strong expectation that a table would have the same projection for every feature - so you can't really use sf for this - there's no sf structure for mixed projections, you have to manage it yourself - usually I'd think that reprojection is fast enough to do it on the fly like this
    – mdsumner
    Sep 24, 2019 at 0:11


Browse other questions tagged or ask your own question.