Search code examples
rimage-processingcropraster

Crop raster with polygon in R: Error extent does not overlap

I want to crop a raster stack using a polygon shapefile i made in ArcGIS, however I get error that extent does not overlap.

First I create the raster stack:

test1 < stack("C:/mydir/test1.tif")

define projection

myCRS <- test1@crs

then read shapefile

myExtent <- readShapePoly("C:/mydir/loc1.shp", verbose=TRUE, proj4string=myCRS)

Crop

myCrop <- crop(test1, myExtent)
Error in .local(x, y, ...) : extents do not overlap

I have searched for a solution, but i only find that projection can be the problem, however they are definetly both in the same CRS:

> test1$test1.1
class       : RasterLayer 
band        : 1  (of  4  bands)
dimensions  : 10980, 10980, 120560400  (nrow, ncol, ncell)
resolution  : 10, 10  (x, y)
extent      : 6e+05, 709800, 5690220, 5800020  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=31 +datum=WGS84 +units=m +no_defs +ellps=WGS84 
+towgs84=0,0,0 
data source : C:\mydir\test1.tif 
names       : test1.1 
values      : 0, 65535  (min, max)

> myExtent
class       : SpatialPolygonsDataFrame 
features    : 1 
extent      : 499386.6, 517068.2, 6840730, 6857271  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=31 +datum=WGS84 +units=m +no_defs +ellps=WGS84 
+towgs84=0,0,0 
variables   : 2
names       :    Shape_Leng,    Shape_Area 
min values  : 67444.6461177, 283926851.657 
max values  : 67444.6461177, 283926851.657
Solution

  • The message is pretty self explanatory. Your extent do not overlap... here how to check it:

    library(raster)
ext.ras <- extent(6e+05, 709800, 5690220, 5800020)
ext.pol <- extent(499386.6, 517068.2, 6840730, 6857271)


plot(ext.ras, xlim = c( 499386.6,709800), ylim= c(5690220,6857271), col="red")
plot(ext.pol, add=T, col="blue")

    I've created extent object from data in your question. You see one extent in the top left corner and the other in the bottom right. Have you tried reading both files in QGIS, you could probably easily see it.

    enter image description here

    If they really are suppose to overlap, than I would suspect the way you read your shapefile. Instead of 

    myExtent <- readShapePoly("C:/mydir/loc1.shp", verbose=TRUE, proj4string=myCRS)

    use :

    library(rgdal)
myExtent <- readOGR("C:/mydir","loc1.shp")
myExtent <- spTRansform(myExtent, CRS(proj4string(test1)))