Clip raster with shapefile and set a value outside the mask

Question

I want to clip a raster using a shapefile with the command:

gdalwarp -cutline INPUT.shp INPUT.tif OUTPUT.tif

From this answer. However I want to set the pixels outside the mask to a particular value (-9999 for example). I want to do it with gdalwarp so I can perform the same process in several files (In other words I don't want to use QGIS, ArcGIS, or any similar package). A Python programmatic approach would also be desirable, as currently I am calling gdalwarp from a Python script.

Answer 1

You need to use the -dstalpha option to gdalwarp e.g.:

gdalwarp -cutline INPUT.shp -crop_to_cutline -dstalpha INPUT.tif OUTPUT.tif

This will add an alpha band to the output tiff which masks out the area falling outside the cutline.

P.S. duplicate question

Answer 2

It is -dstnodata you are asking for, not -dstalpha:

gdalwarp -cutline -dstnodata -9999 INPUT.shp INPUT.tif OUTPUT.tif

gdalwarp sets the area outside the -cutline polygon to -dstnodata, which can be a specified value (or one for each band), the nodata values from the input data set, or, if none of the two is available, some default value, which however is not set as nodata value on the output dataset. I don't know how GDAL determines that default value, but for Byte data, it is zero.

An alternative solution that was proposed here "is to use gdal_translate with -projwin to extract to shape envelope and then gdal_rasterize with -i to burn nodata where shape is missing." This solution also circumvents raster shifting or resampling.

But all of these are using GDAL command line tools.

Answer 3

Today, the pythonic solution would be to use rasterio. For clipping a raster with a polygon, they have an example in their docs.

Answer 4

Set one image to the values you want, copy the full file name of said raster, insert full file name into script below, and drop this script in Python box. Hit enter+enter and sit back and relax.

mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd, "")[0]
refLayer = arcpy.mapping.ListLayers(mxd, "ReplaceWithNameofRaster.tif", df)[0]
for Layer in arcpy.mapping.ListLayers(mxd, "", df):
    arcpy.mapping.UpdateLayer(df, Layer, refLayer)