How to create a GeoTIFF with PyGMT?

Q&A PyGMT Q&A
#1

Hello!

I’m new to GMT and PyGMT, and I’m trying to create a GeoTIFF export of a bathymetry dataset.

I can plot the data just fine as follows:

fig = pygmt.Figure()
fig.coast(region=region, projection="M20c", land="#666666")
fig.grdimage(
    grid=grid.bathymetry,
    cmap="viridis",
    nan_transparent=True,
    # img_out="out.tif",
)
fig.basemap(frame=True)
fig.colorbar(frame='af+l"bathymetric depth [m]"')
fig.show()

image
image2594×5572 484 KB

But when I try to use the out_image parameter to write a GeoTIFF file with

fig = pygmt.Figure()
fig.grdimage(
    grid=grid.bathymetry,
    cmap="viridis",
    nan_transparent=True,
    img_out="out.tif",
)

I get the following error:

grdimage [ERROR]: Option -A: No output argument allowed
grdimage [ERROR]: Option -A: Must provide an output filename for image
---------------------------------------------------------------------------
GMTCLibError                              Traceback (most recent call last)
Cell In[194], line 3
      1 fig = pygmt.Figure()
      2 # fig.coast(region=region, projection="M20c", land="#666666")
----> 3 fig.grdimage(
      4     grid=grid.bathymetry,
      5     cmap="viridis",
      6     nan_transparent=True,
      7     img_out="out.tif",
      8 )
      9 # fig.basemap(frame=True)
     10 # fig.colorbar(frame='af+l"bathymetric depth [m]"')
     11 # fig.show()
     12 # fig.savefig('out.tif')

File /opt/conda/envs/pygeo/lib/python3.11/site-packages/pygmt/helpers/decorators.py:598, in use_alias.<locals>.alias_decorator.<locals>.new_module(*args, **kwargs)
    591     msg = (
    592         "Parameters 'Y' and 'yshift' are deprecated since v0.8.0. "
    593         "and will be removed in v0.12.0. "
    594         "Use Figure.shift_origin(yshift=...) instead."
    595     )
    596     warnings.warn(msg, category=SyntaxWarning, stacklevel=2)
--> 598 return module_func(*args, **kwargs)

File /opt/conda/envs/pygeo/lib/python3.11/site-packages/pygmt/helpers/decorators.py:738, in kwargs_to_strings.<locals>.converter.<locals>.new_module(*args, **kwargs)
    736             kwargs[arg] = separators[fmt].join(f"{item}" for item in value)
    737 # Execute the original function and return its output
--> 738 return module_func(*args, **kwargs)

File /opt/conda/envs/pygeo/lib/python3.11/site-packages/pygmt/src/grdimage.py:193, in grdimage(self, grid, **kwargs)
    190     kwargs["I"] = stack.enter_context(shading_context)
    192 fname = stack.enter_context(file_context)
--> 193 lib.call_module(
    194     module="grdimage", args=build_arg_string(kwargs, infile=fname)
    195 )

File /opt/conda/envs/pygeo/lib/python3.11/site-packages/pygmt/clib/session.py:530, in Session.call_module(self, module, args)
    526 status = c_call_module(
    527     self.session_pointer, module.encode(), mode, args.encode()
    528 )
    529 if status != 0:
--> 530     raise GMTCLibError(
    531         f"Module '{module}' failed with status code {status}:\n{self._error_message}"
    532     )

GMTCLibError: Module 'grdimage' failed with status code 72:
grdimage [ERROR]: Option -A: No output argument allowed
grdimage [ERROR]: Option -A: Must provide an output filename for image

These two errors seem to contradict each other, which is confusing. Does anyone have ideas about how to debug the issue further? Or can anyone suggest a workaround? I don’t mind using other software if it’s not currently possible with pygmt.

Thanks!
Oliver

#2

Hello @OliverEvans96,

Welcome to the GMT forum :slightly_smiling_face:!

I can reproduce this issue. As the related GMT code works, I opened an issue in the PyGMT repository on GitHub, please see: https://github.com/GenericMappingTools/pygmt/issues/2599

I have only little experience with the img_out parameter of Figure.grdimage. Maybe you can use Figure.savefig (https://www.pygmt.org/latest/api/generated/pygmt.Figure.savefig.html) to save the grid as a tif file as a workaround? But I am not sure whether the results are identical and this works for your needs.

import pygmt

fig = pygmt.Figure()
fig.grdimage(grid="@earth_age_01d_g")
fig.savefig(fname="test_pygmt_savefig.tif")
#3

Hello,

I found this discussion when I was looking for information on making a GeoTIFF from PyGMT. In classic GMT, I have long used this command to generate a GeoTIFF file from a EPS file:
gmt psconvert -A -P -Tt -W+g $out

I see there is a pygmt.Figure.psconvert call that wraps the GMT psconvert, but it does not seem to have the -W+g option to make a GeoTIFF.

I tried the `fig.savefig(fname=“out.tif”) call but it made a regular TIFF image without the geolocation information of a GeoTIFF. I tried with the new PyGMT 0.10.0 release that just dropped, but it seems to be the same. Am I missing something, or is the GeoTIFF option not yet supported?

#4

I think it’s not supported yet, most likely because the PyGMT developers don’t use GeoTIFF, but it should be easy to support, so please open an issue report to the PyGMT repository.

1 Like