1

I am trying to crop a raster by a polygon in R to the state of Sergipe-Brazil. My raster is a raster brick criated by 7 .tif files which are EVI (vegetation index) time series. Since the EVI archives are larger then the state I want, I'm trying to crop it by using a Polygon shapefile of Sergipe. Here is the script:


# Loading the packages
library(raster)
library(sp)
library(rgdal)
library(RCurl)
library(bitops)
library(hddtools)
library(XML)
library(zoo)

# Creating the rasterbrick
evi <- list.files(pattern = '^.*\\.tif$')
evistack <- stack(evi)
evibrick <- brick(evistack)

# Converting from -1 to 1 (EVI)
evit1<-evibrick/10000

#Saving the RASTER BRICK
writeRaster(evit1, filename="evi_2000t1.tif", format="GTiff", overwrite=TRUE)

# Importing the shapefile of Sergipe
se <- readOGR(dsn = "C:/ArquivosR/EVI/Dados_anuais/2000", layer = "Estado")

### Transforming the CRS os the shapefile
SE.utm <- SpatialPolygons(list(Polygons(list(Polygon(se)), "1")))
proj4string(SE.utm) = CRS("+proj=utm +zone=24 +south +datum=WGS84 +units=m +no_defs")
SE.utm@bbox; SE.utm@proj4string

### CRS GEO (GeoTiff)
SE.geo<- spTransform(SE.utm, CRS("+proj=longlat +datum=WGS84"))
SE.geo@bbox; SE.geo@proj4string

# Reprojecting the RasterLayer
projection(evit1) <- CRS("+proj=longlat +datum=WGS84")
alt.utm <- projectRaster(evit1, crs="+proj=longlat +datum=WGS84")

### Cropping
evit1.crop <- crop(alt.utm, extent(SE.geo), snap='out')
image(evit1.crop, asp=1)

I've changed the CRS os the shapefile to let both rasterbrick and shapefile in the same coordinate system. BUT, when I run it:

cr <- crop(alt.utm, extent(SE.geo), snap="out")

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

Checking out the extent of the files I got this:

extent(alt.utm)
class       : Extent 
xmin        : -38.24549 
xmax        : -36.39112 
ymin        : -11.56996 
ymax        : -9.51354 

extent(SE.geo)
class       : Extent 
xmin        : -37.44375
xmax        : -37.44375 
ymin        : -10.58074 
ymax        : -10.58074 

I'll aprecciate any help!

  • 1
    SE.geo looks degenerate: its xmin = xmax and ymin = ymax. – J Kelly Apr 3 '16 at 23:09
1

It looks like you have a single point in the shapefile you are trying to crop to and the error message isn't exactly right. Lets' try and simulate that with a reproducible example:

b <- brick(system.file("external/rlogo.grd", package="raster"))
p=data.frame(x=50,y=35) # a point in the middle of b
coordinates(p)=~x+y

Then:

> extent(b)
class       : Extent 
xmin        : 0 
xmax        : 101 
ymin        : 0 
ymax        : 77 
> extent(p)
class       : Extent 
xmin        : 50 
xmax        : 50 
ymin        : 35 
ymax        : 35 
> crop(b, extent(p), snap="out")
Error in .local(x, y, ...) : extents do not overlap

Perhaps the crop algorithm is computing the area of overlap, getting zero, and assuming the areas don't overlap. Anyway, its a minor bug in crop, I'd say. To check, lets just make sure it works with a larger extent:

> p2=data.frame(x=c(50,51),y=c(35,38))
> coordinates(p2)=~x+y
> crop(b, extent(p2), snap="out")
class       : RasterBrick 
dimensions  : 3, 1, 3, 3  (nrow, ncol, ncell, nlayers)
resolution  : 1, 1  (x, y)
extent      : 50, 51, 35, 38  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=merc +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : in memory
names       : red, green, blue 
min values  : 155,   165,  222 
max values  : 157,   168,  224 

Of course we can't be sure this is your situation without your data, but it looks like it. Track back your SE.geo object and see what the original is like:

se <- readOGR(dsn = "C:/ArquivosR/EVI/Dados_anuais/2000", layer = "Estado")

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.