R plant using polygon holding cells along the border
I'm trying to crop this bitmap over Italy, but the output seems to skip some cells along the border. See Areas highlighted in red in the image below:
How can I keep all the cells that cross borders?
Below is my script:
library(raster)
# Load data
x <- raster("x.nc")
IT <- getData(name = "GADM", country = "Italy", level = 0)
# Mask and crop
x_masked <- mask(x, IT)
x_masked_cropped <- crop(x_masked, IT)
# Plot
plot(x_masked_cropped)
plot(IT, add = T)
source to share
You can identify all the raster cells crossing Italy and set the remaining, that is, non-intersecting pixels, to NA. Do not forget to return cells with corresponding weights via cellFromPolygon(..., weights = TRUE)
- otherwise only those cells that have a center within Italy will be returned (see also ?raster::extract
).
## identify cells covering italy and set all remaining pixels to NA
cls <- cellFromPolygon(x, IT, weights = TRUE)[[1]][, "cell"]
x[][-cls] <- NA
plot(trim(x))
plot(IT, add = TRUE)
source to share
Here's one way to do it. We create a raster bitmap layout with gdalUtils::gdal_rasterize
using at=TRUE
to ensure that a value of 1 is burned in all cells affected by the Italy polygon. gdal_rasterize
refers to files on disk, so write IT
to a file with OGR support first.
library(gdalUtils)
library(rgdal)
x_crop <- crop(x, IT)
writeOGR(IT, tempdir(), f <- basename(tempfile()), 'ESRI Shapefile')
gdal_rasterize(sprintf('%s/%s.shp', tempdir(), f),
f2 <- tempfile(fileext='.tif'), at=T,
tr=res(x_crop), te=c(bbox(x_crop)), burn=1,
init=0, a_nodata=0, ot='Byte')
plot(x_crop*raster(f2)) # multiply the raster by 1 or NA
plot(IT, add=TRUE)
source to share