13

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:

  1. I need a way to get a Raster out of the extracted list again or
  2. 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?

3
  • What if it is a 4 bands rgbi raster? Bands are lost so far...
    – Doris
    Jan 25, 2019 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, 2019 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, 2019 at 15:17

3 Answers 3

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)
9
  • 4
    extract() does perform on the fly reprojection, but crop() does not. That might account for some confusion
    – mdsumner
    Apr 8, 2014 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, 2015 at 7: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. Dec 11, 2015 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)) Feb 1, 2018 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. Feb 1, 2018 at 21:01
11

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

mask(raster, poly_shape)

works without errors and returns what I searched for.

8
  • 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. Apr 8, 2014 at 18:08
  • Use projectExtent() to get just the extent for cropping the raster.
    – mdsumner
    Apr 8, 2014 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). Apr 8, 2014 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 9, 2014 at 11:05
  • 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, 2014 at 12:33
-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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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