Hi all,
Condensed problem: I want to export a georeferenced geotiff raster with contours.
I have a bathymetry grid (geotiff). The coordinates are in projected units.
If I use
gmt grdimage -Abat_gmt.tif bat.tif -I+d
I get a nice color-relief map with shading. The coordinates are preserved and I can plot it in e.g. QGIS.
- NaN values are represented with the value 99999 in my grid. How can I tell
grdimageto use this as a NaN-value?
I would like to plot contours as well. That means I need to bake a GMT cake, calling e.g. grdcontour and other modules (as needed). After applying contours I want to export it, keeping the proper coordinates, like the result in the grdimage call above.
Any tips on how doing this?
I tried the psconvert approach (https://docs.generic-mapping-tools.org/latest/psconvert.html#w), adding a %%PROJ line in the .ps-file to specify the projection that should be used. But I’m, confused what to use when the grid is in already in projected units.
To test with geodetic coordinates, I’ve inverse projected the grid, and did a grdimage, and then injected:
%%PROJ utm 6.663782868 8.589722967 73.280813537 73.680863706 238337 295263 8149712 8187663 +proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
in the first line of test.ps and ran gmt psconvert -Tt -W+g test.ps (this creates two files; test.tif and test.tiff - should it do that?)
$ gmt psconvert -Tt -W+g test.ps -V
psconvert [INFORMATION]: Processing test.ps...
psconvert [INFORMATION]: Find HiResBoundingBox ...
psconvert [INFORMATION]: Figure dimensions: Width: 443.772 points [15.6553 cm] Height: 607.74 points [21.4397 cm]
psconvert [INFORMATION]: [79.236 72 523.008 679.74]...
psconvert [INFORMATION]: Convert to TIF...
psconvert [INFORMATION]: width = 2533 height = 1850 X res = 22.473747 Y res = 20.525149
psconvert [INFORMATION]: Wrote world file test.tfw
psconvert [INFORMATION]: Proj4 definition: +proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
Input file size is 2533, 1850
0...10...20...30...40...50...60...70...80...90...100 - done.
psconvert [INFORMATION]: The gdal_translate command:
gdal_translate -mo TIFFTAG_XRESOLUTION=300 -mo TIFFTAG_YRESOLUTION=300 -a_srs '+proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs' -co COMPRESS=LZW -co TILED=YES 'test.tif' 'test.tiff'
The resulting tiff has reasonable y-coordinates, but the x-coordinates are way off:
$ gmt grdinfo test.tif
test.tif: Title: Grid imported via GDAL
test.tif: Command:
test.tif: Remark:
test.tif: Pixel node registration used [Cartesian grid]
test.tif: Grid file format: gd = Import/export through GDAL
test.tif: x_min: -96523725866.8 x_max: -96523668940.8 x_inc: 22.4737465456 name: x n_columns: 2533
test.tif: y_min: 8149712 y_max: 8187683.52515 y_inc: 20.525148729 name: y n_rows: 1850
test.tif: v_min: 24 v_max: 255 name: z
test.tif: scale_factor: 1 add_offset: 0
= Other approaches
== mapnik
I read some blog posts about mapnik; seemed to be a nice solution. You create an .xml-file with all the layers you want (tif’s, shapefiles etc) and mapnik renders a final nice image with all layers combined.
I tried this, but setting up mapnik, for me, was a terrible experience. All kind of python and gdal plugins needing to be manually compiled. It made a mess of my system, with no result.
== gdal
I tried using a workflow based on gdal; 1)make color relief, 2)make hillshade, 3)combine these, 4)extract contours and 5)burn contours into the final raster.
This gives an OK result, but the contours becomes very rough and blocky. This is related to that the finest features cannot be thinner/finer than the raster resolution. One solution is to resample the grid to an artificial higher resolution, so that the contours become nicer, but this cause huge rasters.
gdal has many features that I’ve neved (had to) used, but they are quite handy to have in the toolbox. See below for the gdal workflow I used.
#make cpt
gmt makecpt -Cabyss -T-4000/-2000/50 -N | awk '{print $1, $2}' | sed 's/\// /g' > col.cpt
#resample grid if necessary
#make a color-relief dem based on the cpt created above
gdaldem color-relief -b 1 bat.tif col.cpt bat_col.tif
#make shade
gdaldem hillshade bat.tif bat_shade.tif
#gamma hillshade
gdal_calc.py -A bat_shade.tif --outfile=bat_shade_gamma.tif --calc="uint8(((A / 255.)**(1/0.5)) * 255)" --overwrite
# overlay shade and color-relieff
gdal_calc.py -A bat_shade_gamma.tif -B bat_col.tif --allBands=B \
--calc="uint8( ( \
2 * (A/255.)*(B/255.)*(A<128) + \
( 1 - 2 * (1-(A/255.))*(1-(B/255.)) ) * (A>=128) \
) * 255 )" --outfile=bat_final.tif --overwrite
#extract contours
gmt grdcontour bat.tif -D%c.gmt -C5 -Q20
#add OGR-header. gdal requires this. @Je is the epsd code of the projection
sed -i '1i # @VGMT1.0 @GLINESTRING # @Je25833' C.gmt
#burn in the contours to the raster
gdal_rasterize -b 1 -b 2 -b 3 -burn 0 -burn 0 -burn 0 C.gmt bat_final.tif
#export to tiles, if wanted
#gdal2tiles.py bat_final.tif