Tuesday 26 March 2019

r - Cropping a raster by a polygon - cells missing that are partially outside the polygon


I am trying to crop a raster in RStudio using the crop() function in the Raster package.



I have used the following code:


studyarea<-crop(chlorophyll,mcp2015)

Where 'chlorophyll' is the raster, and 'mcp2015' is a 'SpatialPolygon'


I plotted the output and noticed a problem:


plot(studyarea)
plot(mcp2015, add=TRUE)

enter image description here


The cropped raster is missing cells that are, at least in part, inside the SpatialPolygon. I have used this exact same code, but cropped the raster by a different SpatialPolygon, and I did not have this problem.



I also checked the extent of the cropped raster and the SpatialPolygon, and they are indeed different (ruling out a plotting issue).


>mcp2015

class : SpatialPolygons
features : 1
extent : -130.4485, -129.095, 50.64663, 51.22328 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

studyarea


class : RasterLayer
dimensions : 13, 33, 429 (nrow, ncol, ncell)
resolution : 0.04166667, 0.04166667 (x, y)
extent : -130.4583, -129.0833, 50.66666, 51.20833 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : in memory
names : Chlorophyll.Concentration..OCI.Algorithm
values : 0.444908, 20.26956 (min, max)

Any ideas?




Answer



You need to use snap = 'out', which is one of several clipping options related to the raster cell extent and centre. Also check the crop against the input, you may need to look at AlignExtent() if they don't match up.


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...