Search code examples
rgisrastershapefilemercator

Get mercator coordinates from a list of extent objects in R

I have a shape file with several polygons, and I need to get the mercator coordinates per every polygon in a different file OR in a list file. Here, the list with the extent from each polygon from my shapefile

a<-extent(659302.4, 659802.4, 9860325, 9860825)
b<- extent(647663, 648163, 10033987, 10034487)
c<-extent(609015.1, 609515.1, 9871028, 9871528)
ab<-c(a,b,c)

I need to generate a list of mercator coordinates for each extent in the list OR a file that contains the coordinate per each file. My only option so far has been to rasterize each extent, assign CRS coordinates and then transform to mercator, as follows

bb<-raster(b)
crs(bb)<- "+proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
mbb<-projectRaster(bb, crs="+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0
                   +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs")

mbb class : RasterLayer dimensions : 12, 12, 144 (nrow, ncol, ncell) resolution : 50, 50.3 (x, y) extent : -8869218, -8868618, 34170.81, 34774.41 (xmin, xmax, ymin, ymax) coord. ref. : +proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs

The only problem with this method is that is very slow, specially when I need to get the mercator coordinates of the extent of hundreds of polygons.

EDIT: What I want is to change the projection from UTM to Mercator, I am thinking in a function that allow to do this for the extent list instead of doing individually as showed in the example above.

Solution

  • If you want to change the projection of extent (not raster), you could use:

    library(raster)
library(dplyr)
library(purrr)

a<-extent(659302.4, 659802.4, 9860325, 9860825)
b<- extent(647663, 648163, 10033987, 10034487)
c<-extent(609015.1, 609515.1, 9871028, 9871528)
ab<-c(a,b,c)

ab <- map(ab, function(x) {as(x,'SpatialPolygons')}) %>%
      map(.,function(x){x@proj4string <- CRS("+proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0");
                        spTransform(x,CRS("+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs"))})

# back to extent
new_extent <-  ab %>% map(extent)