0

I want to clip a raster in R based on a polygon. The analog function in arcgis that does the same is Data Management Tools -> Raster -> Raster Processing -> Clip.

What I have until now is:

library(raster)

set.seed(123)

r <- raster()
res(r) <- c( 0.016666, 0.016666 )

extent <- shapefile("colombia_2000.shp")

r <- crop(r, extent, snap='near')
r <- setValues(r, sample(x=c(0,1),size=ncell(r), replace = T))

undata_mask <- mask(r, extent)

You can find the "extent" file in the following link:

https://drive.google.com/open?id=1HLTAkqClWNzVRzWLbTvfp4Esqyc8zMRl

The polygon (extent) I am using represents the area between 1000 an 2000 meters above see level in the Andean Mountain range in Colombia, and as you can see there are some internal patches that are not covered with the polygon, which belongs to the mountain peaks that exceed this altitud range.

When I perform the crop and mask functions in R, those internal patches remain untouched as if they belonged to the polygon. When I perform the calculation in ArcGis, the internal patches are indeed clipped, which is what I want to achieve in R.

How Can I get rid of this internal patches?

Here are some images of What I just said.

The first one is with R and the second one is with ArcGIS

enter image description here

enter image description here

  • You should probably warn people that's an 893Mb raster, which is quite a download. Can you make a much, much smaller example? Like under a megabyte? – Spacedman Nov 7 '17 at 17:05
  • Hi, I am sorry for that extra big raster. I just edit the post and add a reproducible raster instead of the former. – Natalia Arcila Marin Nov 7 '17 at 17:22
1

First, If you want to make a reproducible example... Set resolution before than values (and don't use functions names as object name):

library(raster)

set.seed(123)

r <- raster()
res(r) <- c( 0.016666, 0.016666 )
undata <- setValues(r,values = sample(x=c(0,1),
                    size=ncell(r),replace = T), extent = c( -180, 180, -90, 90))

extnt <- shapefile ("path_to/colombia_2000.shp")
crs.geo  <-CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
crs(extnt) <- crs.geo

I've been suffering the same issue... I tried to use rasterize() in the past, but I got the same problem. Could be a bug o a misinterpretation of neighbours. I don't know the exact reason.

With your data:

enter image description here

I finally decided to take a workaround in QGIS (prior to use only open source software, also, you can make a script to program this with pyqgis).

Create a blank raster using the original one and multiplying it with raster calculator with 0.

"undata@1"*0

After that, in the vector layer, create a new field -decimal number with at least one decimal- with a unique number to all polygons... In this case, 1.1.

enter image description here

Use this field to rasterize the polygon (that's the reason to use a decimal number field type):

enter image description here

Finally, back in R:

undata_crop <- crop(undata, extnt, snap='near')
workaround <- raster("path_to/rasterized.tif")
workaround_crop <- crop(workaround, extnt, snap='near')
undata_mask <- mask(x=undata_crop,mask= workaround_crop ,maskvalue=0)

enter image description here

You can also use undata_crop as input to avoid to crop workaround (but, if undata is actually an input raster, is faster to work directly with it rather than cropping it in R to export to QGIS)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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