1

I have this raster (land-and-water.tif):

enter image description here

And this one (land-only.tif, black = no data):

enter image description here

And I need this (black = no data):

enter image description here

I'm trying to do it with GDAL. Here's what I tried so far:

gdaltindex clipper.shp land-only.tif
gdalwarp -cutline clipper.shp -dstalpha land-and-water.tif water-only.tif

But the clipper.shp is empty, and obviously the resulting water-only.tif is the same as land-and-water.tif . I think I need to tell gdaltindex what is the no-data value in my land-only.tif , but I didn't find how in the doc.

How can I do that?

1

as far as I can estimate that, several processing steps would be neccesary to achieve what you are looking for. In a first step, you should preprocess your land-only raster and throw away any channels and pixel values, to be able to convert it to a simple extents-shapefile without killing your PC in a second step:

gdal_calc.py --calc=1 --outfile=land-only-extents.tif -A land-only.tif  
gdal_polygonize land-only-extents.tif land-only.shp

What you want, is the difference of land-only-extents and the extents of land-and-water.tif. The reason is, if you use gdal_warp with your land only extents like you proposed, you will get your land-only.tif – again. Therefore you also need the extents of your land-and-water.tif. If this is a rectangle, you can use gdaltindex as you already mentioned. If not, use the same process again. Having that, there exists the possibility to calculate the difference of both with ogr2ogr. This source indicates how to do that, and for your case it should look something like that (untested):

ogr2ogr difference.shp land-and-water.shp -dialect SQLite -sql "SELECT ST_Difference(land-and-water.geometry, land-only.geometry) AS geometry FROM land-and-water, 'land-only.shp'.land-only"

Some further information on ogr sql can be found here. This should give you a clipping layer that should work somehow like that:

gdalwarp -cutline difference.shp -dstalpha land-and-water.tif water-only.tif

If you think nodata-values might be a problem, or maybe have to fight with coordinate systems, you can always set that with gdal_translate:

gdal_translate -a_srs <source srid> -s_srs <dest srid> -a_nodata <nodata value> input_raster.tif output_raster.tif

Please let us know if that helped or worked for you.

  • Thanks for the detailed answer! I finally found a solution using Mapnik's compositing and a water shapefile from Natural Earth. – Tim Autin Sep 19 '18 at 10:57

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.