Figure.grdimage: "img_out" ("A") does not work

[yvonnefroehlich]

Related to this post in the GMT Forum: https://forum.generic-mapping-tools.org/t/how-to-create-a-geotiff-with-pygmt/4113

The img_out (A) parameter of Figure.grdimage does not allow to save a grid to a raster format (.bmp, .gif, .jpg, .png, or .tif) as proposed by the documentation https://www.pygmt.org/dev/api/generated/pygmt.Figure.grdimage.html:

import pygmt

fig = pygmt.Figure()
fig.grdimage(grid="@earth_age_01d_g", img_out="test_pygmt_grdimage.tif")

The relatd GMT code works (see https://docs.generic-mapping-tools.org/dev/grdimage.html#a):

gmt begin
    gmt grdimage @earth_age_01d_g -Atest_gmt_grdimage.tif
gmt end

In PyGMT, a figure can be saved as a .bmp, .gif, .jpg, .png, or .tif file via Figure.savefig (https://www.pygmt.org/dev/api/generated/pygmt.Figure.savefig.html). But I am not sure whether this is identical to using img_out of pygmt.Figure.grdimage.

import pygmt

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

Error message

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 

System information

PyGMT information:
  version: v0.9.1.dev98
System information:
  python: 3.11.4 | packaged by conda-forge | (main, Jun 10 2023, 17:59:51) [MSC v.1935 64 bit (AMD64)]
  executable: C:\ProgramData\Anaconda3\envs\pygmt_env_dev\python.exe
  machine: Windows-10-10.0.19045-SP0
Dependency information:
  numpy: 1.24.3
  pandas: 2.0.2
  xarray: 2023.1.1.dev17
  netCDF4: 1.6.2
  packaging: 23.1
  contextily: 1.3.0
  geopandas: 0.13.2
  IPython: 8.14.0
  rioxarray: 0.14.1
  ghostscript: 9.54.0
GMT library information:
  binary version: 6.4.0
  cores: 4
  grid layout: rows
  image layout: 
  library path: C:/ProgramData/Anaconda3/envs/pygmt_env_dev/Library/bin/gmt.dll
  padding: 2
  plugin dir: C:/ProgramData/Anaconda3/envs/pygmt_env_dev/Library/bin/gmt_plugins
  share dir: C:/Program Files (x86)/gmt6/share
  version: 6.4.0

GDAL: gdal 3.7.0

[weiji14]

Thanks @yvonnefroehlich for writing up this detailed bug report! I can reproduce the error locally, and am wondering if this is a case of #896, whereby the 'plotting' function grdimage doesn't output a figure, but to a file instead. Will need to look into the clib part of PyGMT to see if this is an issue on the PyGMT side, or with the GMT C API, cc @seisman.

In the meantime, it might be better to recommend users to use other libraries like rioxarray's to_raster if saving to GeoTIFF. Other image formats like bmp/gif/png might be better handled by OpenCV, scikit-image or others.

[seisman]

I tried the following script:

gmt begin map png
	gmt basemap -Rg -JQ15c -Baf
	gmt grdimage @earth_age_01d_p -Atest_out.tif
gmt end

It produces two files: "map.png" and "test_out.tif". "map.png" only has a frame, and "test_out.tif" has no frame. Looking at the source codes, it seems grdimage -A produces an raster image of the current grid, without writing anything to a PS file.

@PaulWessel Could you please explain what the purpose and use case of the -A option, so that we can decide if we should have this feature on the PyGMT side?

[PaulWessel]

Predate wrappers. @joa-quim added it, but cannot imagine it should be available in modern mode?

[PaulWessel]

The purpose was to use grdimage to directly give a raster image from the grid, with no border etc. But clearly, I think it should be disallowed in modern mode since we cannot tell we are building a PS file. In classic mode I think it is fine and backwards compatible. OK with this, @joa-quim ?

[joa-quim]

Why remove? Let me recall that this option let us create GeoTIFF or GeoPDF images independently of the Ghostscript machinery. And if removed it would have to be removed from the grdimage modern mode docs that is the main source of our documentation (the one that opens by default ).

I’m not sure how it works in modern mode. Maybe it needs to have a full file name to prevent it to be written in the session directory.

[PaulWessel]

I agree, just leave it. However, with -A it skips all the PS setup so perhaps @seisman can tell us what confused him - should we improve the docs to make it very clear that -A basically disable the normal PS flow and only writes a raster file and -B etc are not available?

[seisman]

I'm OK with leaving it as it is.

However, with -A it skips all the PS setup so perhaps @seisman can tell us what confused him - should we improve the docs to make it very clear that -A basically disable the normal PS flow and only writes a raster file and -B etc are not available?

What's the difference between

gmt grdimage @earth_age_01d_g -Atest_gmt_grdimage.tif -Baf -JQ10c

and

gmt begin map tif
  gmt grdimage @earth_age_01d_g -JQ10c
gmt end

As I understand it, they're the same, but grdimage -A writes TIFF directly without producing a PS file.

Is it equivalent to calling gdal_translate directly to convert a netCDF file into a TIFF file?

[joa-quim]

Is it equivalent to calling gdal_translate directly to convert a netCDF file into a TIFF file?

Almost, but the first form, the one involving -A, lets us do reprojection, which gdal_translate does not do. And one can apply shading as well in -A.

[PaulWessel]

Also, we should give error if -A is used with -B.

[seisman]

@PaulWessel In PR #2765, I was working on the grdimage's -A option. For me, the new test related to -A works well but then other tests crash (which should not happen).

When I read the grdimage source code, I see that the GMT_IMAGE *Out object (https://github.com/GenericMappingTools/gmt/blob/a8f05c9bc6aed955be8f974928b0abf008d8f42c/src/grdimage.c#L1275) is the one that holds the output image, but it seems it's never freed (I don't see a GMT_Destroy_Data call). Is it intentional or is it a potential bug?

[PaulWessel]

Don't think so. Lots of places were we dont call GMT_Destroy_Data in the module. They are all handled by the garbage collector gmtlib_garbage_collection when the module ends.

[seisman]

For wrappers, the -A option doesn't work for formats that need a driver (e.g., pdf=PDF).

The following bash script works:

gmt begin map
gmt grdimage @earth_relief_01d_g -Cearth -JW0/6i -Atmp.pdf=PDF
gmt image logo.png -Dx0/0+w2c -F+pthin,blue
gmt end show

but the equivalent Python script doesn't:

import pygmt
fig = pygmt.Figure()
fig.grdimage("@earth_relief_01d_g", cmap="earth", projection="W0/6i", img_out="tmp.pdf=PDF")
fig.image("logo.png", position="x0/0+w2c", box="+pthin,blue")
fig.savefig("map.png")

Using img_out="tmp.jpg" or img_out="tmp.png" works.

[joa-quim]

PDF is a vector format, not an image.

[seisman]

PDF is a vector format, not an image.

I have no idea about what it really does, but the docs says:

For other output formats you must append the required GDAL driver. The driver is the driver code name used by GDAL; see your GDAL installation’s documentation for available drivers. Append a +coptions string where options is a list of one or more concatenated number of GDAL -co options. For example, to write a GeoPDF with the TerraGo format use =PDF+cGEO_ENCODING=OGC_BP.

I think it means =PDF is allowed?

[joa-quim]

Funny that it even works in some cases in PyGMT because in Julia

grdimage [ERROR]: Option -A: No output argument allowed

and that comes from the condition

if (API->external) {	/* External interface only */
	if ((n = strlen (opt->arg)) > 0) {
		GMT_Report (API, GMT_MSG_ERROR, "Option -A: No output argument allowed\n");

so I don't see how it can work at all for PyGMT. Is it not detected as an external at that point?

[seisman]

so I don't see how it can work at all for PyGMT. Is it not detected as an external at that point?

In PyGMT, when -Aimg_out.png is given, we call the grdimage module using the following syntax:

gmt grdimage ingrid.nc -A > img_out.png

[PaulWessel]

@joa-quim would know but grdimage -A ends up calling GDAL for writing and I dont know if it is OK there to write a binary image to stdout.

[joa-quim]

I don't remember any details but we certainly pass a file name to GDAL for it to save the image. No idea what redirecting to stdout does.

Tried and it errors on CLI

C:\v>grdimage @earth_relief_01d_g -Cearth -JW0/6i -A > lixo.png
grdimage [ERROR]: Option -A: No output name provided

[seisman]

I don't remember any details but we certainly pass a file name to GDAL for it to save the image. No idea what redirecting to stdout does.

Tried and it errors on CLI

C:\v>grdimage @earth_relief_01d_g -Cearth -JW0/6i -A > lixo.png
grdimage [ERROR]: Option -A: No output name provided

This doesn't work for GMT CLI. For external wrappers, there are some special handling with > lixo.png. Related codes are:
https://github.com/GenericMappingTools/gmt/blob/29001d863941e6e4fd3f98a347dba93d63b86250/src/grdimage.c#L257C9-L257C9

https://github.com/GenericMappingTools/gmt/blob/29001d863941e6e4fd3f98a347dba93d63b86250/src/grdimage.c#L1299

[joa-quim]

OK, but actually this complained but worked for me.

julia> gmt("grdimage @earth_relief_01d_g -Cearth -JW0/6i -A > lixo.pdf:PDF")
ERROR 1: Point outside of projection domain
ERROR 1: Point outside of projection domain
ERROR 1: Point outside of projection domain
ERROR 1: Point outside of projection domain

[seisman]

OK, but actually this complained but worked for me.

julia> gmt("grdimage @earth_relief_01d_g -Cearth -JW0/6i -A > lixo.pdf:PDF")
ERROR 1: Point outside of projection domain
ERROR 1: Point outside of projection domain
ERROR 1: Point outside of projection domain
ERROR 1: Point outside of projection domain

Yes, it also works in PyGMT, but it crashes if I call gmt image after this command. In other words, this one complains but works:

import pygmt
fig = pygmt.Figure()
fig.grdimage("@earth_relief_01d_g", cmap="earth", projection="W0/6i", img_out="tmp.pdf=PDF")

This one also works:

import pygmt
fig = pygmt.Figure()
fig.image("logo.png", position="x0/0+w2c", box="+pthin,blue")

but this one crashes:

import pygmt
fig = pygmt.Figure()
fig.grdimage("@earth_relief_01d_g", cmap="earth", projection="W0/6i", img_out="tmp.pdf=PDF")
fig.image("logo.png", position="x0/0+w2c", box="+pthin,blue")

I think the equivalent Julia code is:

> gmt("grdimage @earth_relief_01d_g -Cearth -JW0/6i -A > lixo.pdf:PDF")
> gmt("image logo.png -Dx0/0+w2c -F+pthin,blue")

[joa-quim]

So it crashes when you call image after grdimage -A?
Can´t reproduce it. This works.

julia> gmt("grdimage @earth_relief_01d_g -Cearth -JW0/6i -A > lixo.pdf:PDF")
julia> gmt("psimage lixo.png -Dx0/0+w2c -F+pthin,blue > lixo2.ps")

But this is puzzling me

julia> image("lixo.png", D="x0/0+w2c", F="+pthin,blue", show=1, Vd=1)
        psimage lixo.png  -JX15c/0 -Baf -BWSen -F+pthin,blue -Dx0/0+w2c -P -K > c:\TMP/GMTjl_j.ps
psimage [ERROR]: Option -D option: Modifier +p unrecognized

There is no +p in -D option in the psimage call as can be seen in the generated command.

[joa-quim]

There must be a bug in psimage

C:\v>psimage lixo.png  -R0/1/0/1 -JX7 -Baf -BWSen -F+pthin,blue -Dx0.5c/0.5c+jBL+w6c -P > lixo.ps
psimage [ERROR]: Option -D option: Modifier +p unrecognized

Probably -F is not being parsed

[seisman]

Closing the issue because we won't implement it in PyGMT.