1

I need crop a raster by a grid shapefile to obtain raster files cutted according to grid cells (which are named "PageNumber"). My function is:

g=shapefile("50x50km_final.shp")

f=raster("declivcerrado.tif")

cell=nrow(g)

i=1

for (i in 1:cell){ print(paste("Processando paisagem",i,"de",cell)) tmp=subset(g,g$PageNumber==i) tmp1=crop(f,tmp) writeRaster(tmp1,paste("paisag",i,sep=""),format='HFA',overwrite=T) }

But I got this error:

Error in .local(x, y, ...) : extents do not overlap

The files extent and projection are:

extent(g) class : Extent xmin : -654873.4 xmax : 1381884 ymin : -1416856 ymax : 1056824

extent(f) class : Extent xmin : -60.40702 xmax : -40.48814 ymin : -24.70977 ymax : -2.266753

projection(g) [1] "+proj=aea +lat_1=-2 +lat_2=-22 +lat_0=-12 +lon_0=-54 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs"

projection(f) [1] "+proj=longlat +ellps=GRS80 +no_defs

How can I solve this problem?

  • Start with library(rgdal); g <- spTransform(g, projection(f)); to put the vector data into the same projection as the raster. – mdsumner Jan 17 '17 at 15:41
  • Thank you a lot! I've tryied the transformation and it really works. But would be more interesting for me the opposite, put the raster into the same projection as vector. Is it possible? How? – Thallita Grande Jan 17 '17 at 19:35
  • projectRaster, but this is of course an expensive and destructive resampling process with many complicatedoptions for optimization. ( Try leaflet::addRaster for an interactive view, and just for fun! (though only in Mercator project atm). – mdsumner Jan 17 '17 at 20:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.