[GIS] Extract Raster from Raster using Polygon shapefile in R

linerrastershapefile

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?

Best Answer

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)
Related Question