2

Background: I have a shapefile with multiple polygons over a large area and then I also have GEOTiff image files. I want to crop the TIFF files with each of the polygons. The problem is that when I use raster::intersect to try and do it, I get the error:

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

This is my (pseudo)code for loading the shapefile and selecting a polygon from that shapefile:

 Mineralrights.shape<-readOGR(dsn = "C:/Users/blablabla/Shapefiles",layer = "Mineral Rights")  
 Mineralrights.shape <- spTransform(Mineralrights.shape, "+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs")  
 tempshape <- Mineralrights.shape[Mineralrights.shape$GID == "2097674",]

And this is my (pseudo)code for loading the GEOTiff as a raster:

  Aerialfiletocheck <- raster(paste0(getwd(),"/","2923AD_06_2014_783.tif"))  
  projection(Aerialfiletocheck) <- projection(tempshape)

I don't think the projections is the issue here, since all the source files (shapefile and TIFFS) uses a standard WGS84 ellipse, I just set the projections explicitly for safety. Then call:

  clip <- raster::intersect(Aerialfiletocheck,tempshape)

And this is when I get the error about the extents that do not overlap. To be sure, I tried to simplify the GEOTiff to just its extent and then plotting that over the specific polygon I have selected:

  raam <- projectExtent(Aerialfiletocheck,"+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs")
  raam.shape <- as(raam@extent, 'SpatialPolygons')
  plot(raam.shape)
  plot(tempshape,add=T)

Yielding this plot, where the large rectangle is the extent of the GEOTiff and the funny triangle being the polygon I selected. enter image description here

These are the parameters requested

 > Aerialfiletocheck@extent
    class       : Extent 
    xmin        : 24127.25 
    xmax        : 29297.25 
    ymin        : -3248250 
    ymax        : -3242397 
 > tempshape@bbox
            min       max
    x  23.25377  23.25875
    y -29.33836 -29.33380

Problem: To me it looks like their extents very much overlaps, hence my question: Does the inputs' extent have to be identical for raster::intersect to crop a RGB GEOTiff? That would make no sense to me, because the whole idea of raster::intersect is to "crop" the one with the other, which would imply that their extents will overlap, but not necessarily be identical.

What else have I tried: 1) I have replaced the raster with another polygon and then get the expected result.

2) I have tried base::cropwhich give the same error, which I expected because I think when one input is a raster, intersect just calls crop

3) I have tried rgeos:gIntersection which results in an error about the raster not having a proj4string slot, which makes no sense to me.

4) I have read these posts (and many others): Extracting intersection areas in R

Crop a raster file in R

Extract Raster from Raster using Polygon shapefile in R

GDAL: Convert specific RGB (GeoTIFF) color to transparent

  • Please, add to the post Aerialfiletocheck@extent and tempshape@bbox – aldo_tapia Dec 12 '17 at 11:30
  • @aldo_tapia I see what you did there. It seems the projection of the tiff is in Easting/Westings or polar coordinates or something. Is projection(Aerialfiletocheck) <- projection(tempshape) somehow not working? – Waldo vdM Dec 12 '17 at 12:00
  • Yes, check my answer – aldo_tapia Dec 12 '17 at 12:22
1

You need to use projectRaster(Aerialfiletocheck,tempshape) to change raster's projection... BUT! do not change raster projection, change vector.

If you pretend to project raster, you have two methods to resample information (nearest neighbor and bilinear interpolation). As a resample, new information will be created, so you'll change the original data of your raster.

Instead of this, use spTransform(tempshape, Aerialfiletocheck@crs) to project vector. With the new vector object you can crop, mask, and so on.

  • you beauty! It works. – Waldo vdM Dec 13 '17 at 10:03
  • Can you explain why my initial plot showed that their extents do overlap even though they were in different coordinate systems? – Waldo vdM Dec 13 '17 at 10:05
  • @WaldovdM as function documentation says: projectExtent returns a RasterLayer with a projected extent, but without any values. The CRS used is the same one that shape layer, so the plot shows overlapping features. – aldo_tapia Dec 13 '17 at 10:23

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.