Extract Raster from Raster using Polygon shapefile in R

Question

I am new to R and using the raster package. I have a problem extracting polygons from an existing raster file. If I use

extract(raster, poly_shape)

function on the raster it always creates a list with the data. What I really want is to extract another raster file that I can load with ArcGIS again. After reading a bit more I think the crop function is what I really need. But when I try to use this function

crop(raster, poly_shape)

I get this Error:

Error in .local(x, y, ...) : extents do not overlap
In addition: Warning message:
In intersect(extent(x), extent(y)) : Objects do not overlap

The files raster and poly_shape are the same for both functions. Can you tell me what could be wrong here? Is it even right that the crop function creates another raster and not a list?

EDIT: The extent() function does not work for me. I still get the same error. But I am sure the 2 datasets overlap! With the

extract(raster, poly_shape)

I get the right data from it. Just as a list and not as a raster like I want to have it. I just loaded the datasets in ArcGIS before and they fit very well so I did not check the projection. Now I tried

projection(raster) # "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs"
projection(poly_shape) # "+proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs"

and you can see that the projections do not fit. The extract function seems to be able to automatically transform the files in the right way. I know that because I did the following:

I cut out the exact part of the polygon I extracted in R also in ArcGIS
I calculated the sum of all values of the extracted R polygon (list)
I calculated the sum of all raster cells that I cut out in ArcGIS

The 2 have the exact same result so I guess the conclusion should be that the extract function did work correct. Now I have 2 options I guess:

I need a way to get a Raster out of the extracted list again or
The 2 datasets (raster + poly_shape) need to use the same prjection and the crop function should work

What would you suggest to do here?

What if it is a 4 bands rgbi raster? Bands are lost so far... — Doris, Jan 25 '19 at 14:48
Welcome to GIS SE! As a new user, please be sure to take the short [tour]. Then consider editing your answer to provide additional information and references. See How do I write a good answer? for more info. — Andy, Jan 25 '19 at 15:12
If you have a new question, please ask it by clicking the Ask Question button. Include a link to this question if it helps provide context. - From Review — Erik, Jan 25 '19 at 15:17

score 22 · Answer 1 · answered Apr 07 '14 at 16:59

22

The extract function is behaving exactly as it should. You can force the crop function to use the extent of the polygon and then mask the object to return the exact raster representing the polygon area. If you continue to receive the error it means that your data, in fact, does not overlap.

Please keep in mind that R does not perform "on the fly" projection so, check your projections. You can check if your extents overlap using the "extent()" function.

Here is an example of cropping using the polygon extent then masking the resulting raster using the "rasterized" polygon.

# Add required packages
require(raster)
require(rgeos)
require(sp)

# Create some data using meuse 
data(meuse)
  coordinates(meuse) <- ~x+y
    proj4string(meuse) <- CRS("+init=epsg:28992")
data(meuse.grid)
  coordinates(meuse.grid) = ~x+y
    proj4string(meuse.grid) <- CRS("+init=epsg:28992")
      gridded(meuse.grid) = TRUE    
        r <- raster(meuse.grid) 
          r[] <- runif(ncell(r))

# Create a polygon
f <- gBuffer(meuse[10,], byid=FALSE, id=NULL, width=250, 
                         joinStyle="ROUND", quadsegs=10)   

# Plot full raster and polygon                       
plot(r)
  plot(f,add=T)

# Crop using extent, rasterize polygon and finally, create poly-raster
#          **** This is the code that you are after ****  
cr <- crop(r, extent(f), snap="out")                    
  fr <- rasterize(f, cr)   
    lr <- mask(x=cr, mask=fr)

# Plot results
plot(lr)
  plot(f,add=T)

answered Apr 07 '14 at 16:59

Jeffrey Evans

31,750
2
47
94

4

extract() does perform on the fly reprojection, but crop() does not. That might account for some confusion – mdsumner Apr 08 '14 at 19:12
1

@Jefferey crop() and mask() only clip the raster according to the rectangular extents of the polygon it does not clip it from within the boundary of the polygon. Any idea what commands could clip the raster within the boundary of the given polygon? – csheth Dec 11 '15 at 07:09
1

@Chintan Sheth, for mask to subset within the polygon you need to have a raster representing the values within the polygon. This is why you rasterize the subset polygon and then mask to it. The crop step is to reduce the extent of the raster to the same as the subset polygon, which in the past, if skipped, would throw an extent mismatch error. – Jeffrey Evans Dec 11 '15 at 16:00
spTransform from the sp package (that is sometimes automatically loaded with other spatial R packages) allows to reprojection so that both files are in the same projection eg. good_poly=spTransform(spolygon, CRSobj=crs(raster_file)) – user3386170 Feb 01 '18 at 20:32
@user3386170, Huh? Not sure what you are getting at. This question occurred at a time when the raster package just added "on the fly projection" within some if its functions. These functions had previously thrown an error when projections did not match, but this post was from 2014. You should also be aware of always loading rgdal when using sp::spTransform() as, it adds additional, important, functionality behind the scenes. – Jeffrey Evans Feb 01 '18 at 21:01
My turn for "Huh?". Are you saying that matching projections doesn't matter anymore with recent raster package updates? I just tried to crop a raster with a polygon that had different projections and I ran into this error: Error in .local(x, y, ...) : extents do not overlap. When I transformed the polygon, the crop was performed. I am "getting at" proposing a solution for someone who runs into rasters and polygons with different projections causing problems. – user3386170 Feb 01 '18 at 21:07
Only certain functions eg., extract(). It makes sense that something like crop and mask would not provide on-the-fly projection. I am notable weary of on-the-fly projections and would rather have the error. It keeps people honest. – Jeffrey Evans Feb 01 '18 at 21:21
Still not sure why you don't think it is relevant to guide users to how to resolve issues when polygon and raster are in different projections (if you say on-the-fly can be problematic, I trust your judgement on that). This post is still useful/relevant three years later, e.g. your example/instructions helped me perform this operation on my dataset. – user3386170 Feb 02 '18 at 15:12
The extent of the cropping features had to be within the extent of x for me. It errored with Error in methods::validObject(y) : invalid class “Extent” object: TRUE Seems like a bug. help(crop) reads: "Areas included in y but outside the extent of x are ignored." This doesn't suggest an error would be thrown. – Wassadamo Feb 06 '19 at 09:57

Lars · Accepted Answer · 2016-10-25T08:01:07.733

11

What I actually searched for was the mask() function.

mask(raster, poly_shape)

works without errors and returns what I searched for.

edited Oct 25 '16 at 08:01

answered Apr 08 '14 at 14:52

Lars

2,197
2
18
33

2

Reproject your data so it is in the same projection space. Even in ArcGIS, where on the fly projection is automatic, it is very bad practice to conduct analysis in different projections. With data in different projections you data will not share common extents, which is the error you are receiving. – Jeffrey Evans Apr 08 '14 at 18:08
Use projectExtent() to get just the extent for cropping the raster. – mdsumner Apr 08 '14 at 19:15
To fit the Q&A format of the site this should be placed into the main body of the question as an edit/update (and then add a comment to the answer this is in "reply" to so they know there's more to look at). – matt wilkie Apr 08 '14 at 20:23
@mattwilkie Sorry for not fitting the format but my text was too long to post it as a comment here. @JeffreyEvans I tried the following: projection(raster) = projection(poly_shape) and the other way around projection(poly_shape) = projection(raster) but both ways produce the same error: Error in .local(x, y, ...) : extents do not overlap In addition: Warning message: In intersect(extent(x), extent(y)) : Objects do not overlap. Is there a way how I can see which projection is used on the fly by using the extract() function (because that one obviously works)? – Lars Apr 09 '14 at 11:05
I solved this by exporting the 2 datasets from ArcMap again using the current dataframe as reference. Thank you for your help! – Lars Apr 09 '14 at 20:09
The crop function still is not what I searched for because it extracts the bbox and not just data inside the polygon. Any ideas on how to solve that? – Lars Apr 09 '14 at 20:43
1

What I actually searched for was the mask() function. mask(raster, poly_shape) works without errors and returns what I searched for. – Lars Apr 10 '14 at 12:33
@user3447301 The mask answer seems correct and needs to be more prominent. – geotheory Feb 29 '16 at 09:46

score -1 · Answer 3 · answered Sep 13 '15 at 14:52

-1

Extent works just fine ... I think your extent's Xmin, Xmax, Ymin, and Ymax are different from your raster's X and Y - i.e. they are set opposite

answered Sep 13 '15 at 14:52

ToNoY

573
1
5
14

Extract Raster from Raster using Polygon shapefile in R

3 Answers3

Linked