5

I have a raster of North America and I'd like to remove the larger lakes from it, replacing them with the same no-data value as the original raster.

So I downloaded Natural Earth's 10 meter lakes vector data, which has lakes conveniently ranked by scale rank: 0 for large lakes like the Great Lakes, and larger scale ranks for more minor lakes.

The closest I've come with gdalwarp is:

$ gdalwarp -cutline Data/ne_10m_lakes/ne_10m_lakes.shp \
  -csql "SELECT * FROM ne_10m_lakes WHERE scalerank=0" \
  INPUT.tif lakes-only.tif

This actually does the opposite of what I want, it keeps only the Great Lakes portions of my raster, replacing all other non-lake pixels with the no-data value.

My question: how can I invert the cutline, so instead of getting the area inside the lake boundary, I get whatever is outside it? Can I invert the original Natural Earth .shp file, perhaps?

References

In How to invert polygons (or other options for exterior styling)?, GRASS' r.mask was suggested, which has an inverse mode. Is there something similar with OGR/GDAL?

In How to get an ocean shapefile from a very detailed land area shapefile?, the .shp file was converted to a raster, then inverted, and used as a mask. I'd prefer not to generate a raster of non-lakes because my original TIFF raster is very large (high resolution) and I'd prefer not to generate a raster of the same size. As a stopgap though, I was able to do this:

# Create a not-lakes raster of the same spacing and extent as INPUT
$ gdal_rasterize -i -burn 255 -a_nodata 0 -ot Byte \
  -te `gdalinfo INPUT.tif | egrep "Lower Left|Upper Right"|cut -d"(" -f2 | cut -d")" -f1 | sed 's/,//'| tr "\n" " "` \
  -sql "SELECT * FROM ne_10m_lakes WHERE scalerank=0" \
  -tr `gdalinfo INPUT.tif |grep "Pixel Size"|cut -d"(" -f2|sed 's/[,)]/ /g'` \
  Data/ne_10m_lakes/ne_10m_lakes.shp lakes.tif
# Apply the lakes raster to input as some kind of mask
$ gdal_calc.py -A INPUT.tif -B lakes.tif --outfile OUTPUT.tif --calc='A' \
  --NoDataValue=0

Can something similar be done without creating a gigantic non-lake raster?

4
  • 2
    You can use gdal_rasterize and burn lake polygons into the image with fixed color.
    – user30184
    Commented Dec 10, 2015 at 6:22
  • That would be really neat, using gdal_rasterize to modify the raster without the creation of a secondary lakes-related raster. But the documentation doesn't seem to specify any way to manipulate a raster in this way with gdal_rasterize? It seems to accept one data source and one destination file? Commented Dec 10, 2015 at 6:35
  • @user30184 got it! I thought the "destination filename" was strictly for outputs, but I just tried to give it the input file and what you said worked perfectly! Can you post an answer so I can accept it 😄? Commented Dec 10, 2015 at 6:40
  • @user30184 Write this in an answer: “gdal_rasterize can burn a fixed color for polygons into an existing raster: gdal_rasterize -burn 0 -sql "SELECT * FROM ne_10m_lakes WHERE scalerank=0" Data/ne_10m_lakes/ne_10m_lakes.shp INPUT.tif” so I can accept it! Commented Dec 21, 2015 at 15:42

1 Answer 1

4

gdal_rasterize http://www.gdal.org/gdal_rasterize.html can burn a fixed color for polygons into an existing raster:

gdal_rasterize -burn 0 -sql "SELECT * FROM ne_10m_lakes WHERE scalerank=0" Data/ne_10m_lakes/ne_10m_lakes.shp DESTINATION.tif

According to the documentation the output raster must support update mode access. Which formats support update mode can be checked with gdalinfo by looking at + sign in the driver capabilities. However, updating compressed images may fail or it can lead to sub-optimal result. Uncompressed TIFF is the safest choice but in theory other formats from this list should work as well.

gdalinfo --formats | find "+"

FITS -raster- (rw+): Flexible Image Transport System
HDF4Image -raster- (rw+): HDF4 Dataset
KEA -raster- (rw+): KEA Image Format (.kea)
netCDF -raster- (rw+s): Network Common Data Format
VRT -raster- (rw+v): Virtual Raster
GTiff -raster- (rw+vs): GeoTIFF
NITF -raster- (rw+vs): National Imagery Transmission Format
HFA -raster- (rw+v): Erdas Imagine Images (.img)
ELAS -raster- (rw+v): ELAS
MEM -raster- (rw+): In Memory Raster
BMP -raster- (rw+v): MS Windows Device Independent Bitmap
PCIDSK -raster,vector- (rw+v): PCIDSK Database File
PCRaster -raster- (rw+): PCRaster Raster File
ILWIS -raster- (rw+v): ILWIS Raster Map
SGI -raster- (rw+): SGI Image File Format 1.0
Leveller -raster- (rw+): Leveller heightfield
Terragen -raster- (rw+): Terragen heightfield
ISIS2 -raster- (rw+v): USGS Astrogeology ISIS cube (Version 2)
ERS -raster- (rw+v): ERMapper .ers Labelled
RMF -raster- (rw+v): Raster Matrix Format
RST -raster- (rw+v): Idrisi Raster A.1
INGR -raster- (rw+v): Intergraph Raster
GSBG -raster- (rw+v): Golden Software Binary Grid (.grd)
GS7BG -raster- (rw+v): Golden Software 7 Binary Grid (.grd)
PNM -raster- (rw+v): Portable Pixmap Format (netpbm)
ENVI -raster- (rw+v): ENVI .hdr Labelled
EHdr -raster- (rw+v): ESRI .hdr Labelled
PAux -raster- (rw+): PCI .aux Labelled
MFF -raster- (rw+v): Vexcel MFF Raster
MFF2 -raster- (rw+): Vexcel MFF2 (HKV) Raster
BT -raster- (rw+v): VTP .bt (Binary Terrain) 1.3 Format
LAN -raster- (rw+v): Erdas .LAN/.GIS
IDA -raster- (rw+v): Image Data and Analysis
GTX -raster- (rw+v): NOAA Vertical Datum .GTX
NTv2 -raster- (rw+vs): NTv2 Datum Grid Shift
CTable2 -raster- (rw+v): CTable2 Datum Grid Shift
KRO -raster- (rw+v): KOLOR Raw
ROI_PAC -raster- (rw+v): ROI_PAC raster
ISCE -raster- (rw+v): ISCE raster
ADRG -raster- (rw+vs): ARC Digitized Raster Graphics
SAGA -raster- (rw+v): SAGA GIS Binary Grid (.sdat)
PDF -raster,vector- (rw+vs): Geospatial PDF
GPKG -raster,vector- (rw+vs): GeoPackage

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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