1

I am trying to "mask" raster cells that are outside of a polygon layer, setting them to "NoData". Also, I need to use GDAL, rather than R which I normally use, because the raster is too large to fit in memory.

So far, my workflow is -

## Burn -9999 into areas outside of polygons
gdal_rasterize -i -b 1 -b 2 -burn -9999 -l pol pol.shp modis_south_copy.tif

## Set -9999 to NoData
gdal_translate -of GTiff -a_nodata -9999 modis_south_copy.tif modis_south_copy1.tif

The first expression sets all cells that are outside of polygon pol.shp in raster modis_south.copy.tif to -9999. The second expression creates a new raster modis_south_copy1.tif where the value -9999 is replaced with NoData.

Here is a graphical illustration of the process:

I start with a raster named modis_south_copy.tif, here is an image of the first four bands -

After running the first expression with gdal_rasterize, areas outside pol.shp are replaced with -9999 -

Finally, after the second expression with gdal_translate, the copy modis_south_copy1.tif has -9999 replaced with to NoData -

The input files modis_south.tif and pol.shp are available here - https://www.dropbox.com/sh/b9wtqgpylohuerx/AACyzhzfeGqgthGxMe1ZfW1pa?dl=0

Can I somehow specify that I want to -burn NoData instead of -burn -9999, in gdal_rasterize and, if so, what is the right syntax? (-burn NoData is just a guess and doesn't work)

0

If you like R then I recommend gdalUtils. It interfaces R directly to GDAL, allowing you to use it in a more R-like way; you'll need to write data to disk for GDAL to work with it. The following will mask data to an input geom.

library(gdalUtils)
library(raster)
library(rgdal)

## r raster
## geom SpatialPolygon
## outfile path to save masked data
ezmask <- function(r, geom, outfile, dstnodata='None', ...) {
    tmpd <- tempdir(check=TRUE)
    tmpf <- tempfile(fileext='.tif')
    rgdal::writeOGR(geom, tmpd, 'cutline', driver='ESRI Shapefile', overwrite_layer=TRUE)
    raster::writeRaster(r, tmpf)
    gdalUtils::gdalwarp(tmpf, outfile, cutline=tmpd, cl='cutline', dstnodata=dstnodata, ...)
}

## assuming data is available in home directory
rModis <- raster('~/modis_south.tif')
spPol <- readOGR('~/pol.shp')
ezmask(rModis, spPol, 'masked_modis_south.tif')
rModisMasked <- raster('masked_modis_south.tif')
plot(rModisMasked, ylim=c(34e5, 355e4))

## pass arguments expected by gdalwarp to alter behavior
ezmask(rModis, spPol, 'masked_modis_south.tif', dstnodata=-9999, overwrite=TRUE)

## happy mapping!

modis masked

  • Thanks! this doesn't work for me yet - will see if I can figure out why. The error message is ERROR 5: -tap option cannot be used without using -tr.. – Michael Dorman May 15 at 6:28
  • oh just remove tap. The tr option is to resample to a target resolution--if you need to resample multiple rasters then use tap to ensure they align. I've edited my answer to remove it. – 0mn1 May 16 at 18:05
  • This writes a raster 50 times larger in file size, with all cells either 0 or 'NoData', and at a different CRS and resolution... Not what I needed, but thanks anyway for your time looking into it! – Michael Dorman May 19 at 8:22
  • that's odd. I'll fetch your data today and test it. did you give it a tr argument? that would change resolution, otherwise I don't see how that's possible. the CRS changed because I set it to default to 4326--a bad assumption that I'll remove. – 0mn1 May 21 at 16:18
  • please copy the new function and try it out. I've added a workflow demonstration and plot of results. pass named arguments to alter behavior--you might have inadvertently done this if you passed an argument named tr, which sets a new resolution. – 0mn1 May 21 at 17:08

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.