Setting NaN value in GeoTIFF output with PyGMT

Q&A PyGMT Q&A
1

I am outputting GeoTIFFs from irregularly spaced data (i.e. with big gaps that should be transparent ultimately) as follows:

fig3 = pygmt.Figure()
pygmt.makecpt(cmap='viridis',series=[-15,15])
fig3.grdimage(grid=gridlatlon,cmap=True,nan_transparent=True,
    projection=projection,region=georegion)
fig3.savefig('../../Day1tvg.tiff',crop=True,dpi=300)

This produces a fine output, except the NaN values are all set to 255. Unfortunately, this means when I pull the GeoTIFF into QGIS for example and set the Null value to 255, the maximum valued data points also become transparent. Is there a way to set what the null value is in the geotiff, or perhaps the output data range from 0-254? Thank you

2

try adding transparent=True:

fig3.savefig('../../Day1tvg.tiff',crop=True,dpi=300, transparent=True)
3

That only works for PNG and throws an error with tiffs.

4

Right.

Then you probably need to use grdimage’s -A option. -A has not been implemented as far as I can understand, so if you have to stick to pyGMT you’d have to call grdimage with -A via gmt C api, see below. It may be hard or not, depending on whether you are familiar with command-line GMT. The code below is incomplete, messy and not tested but I did a few things like that.

A separate problem is that the grid data has to be provided as a file in some way, see below.

import pygmt
from pygmt.clib import Session
# no pygmt figure needs to be started to generate cmap and the tiff
cmap_file = 'viridis_based_cmap.cpt' 
pygmt.makecpt(cmap='viridis',series=[-15,15], output=cmap_file)
# the grid data needs to be provided as a file that command-line GMT is able to use. One possible way is to save these to a grid_file.
grid_filename = ... 
# as grid files usually contain a region and a projection, specifying these separately might not be needed.
region = ... # if needed - the region string as W/E/S/N for use with -R, the way command-line GMT understands it
projection = ... # if needed - the projection string for use with -J, the way command-line GMT understands it
tiff_output_file = '../../Day1tvg.tiff'
with Session() as s:
    # -R{region} -J{projection} may be omitted if grid_filename already contains it   
    s.call_module('grdimage', args=f'{grid_filename} -R{region} -J{projection} -E300 -A{tiff_output_file} -C{cmap_file} -Q')

Another way to feed the grid data is using pyGMT internal virtual file mechanism (to avoid explicitly saving it as a grid_file above):

with Session() as s:
    with s.virtualfile_from_grid(gridlatlon) as fin: # gridlatlon as in the example code in the topic
        s.call_module('grdimage', args=f'{fin.name} -E300 -A{tiff_output_file} -C{cmap_file} -Q')

links pygmt.clib.Session.virtualfile_from_grid — PyGMT
grdimage — GMT 6.5.0 documentation

1 Like