Strange behaviour when using Session.call_module().

[timrlaw]

Description of the problem

Hi all. I'm having issues related to sessions, and the order in which seemingly unrelated commands are issued affecting the final plot in strange ways. I think this is related to #1242 and #733.

Essentially what I'm trying to do is set the dimensions of an inset figure to match the dimensions of the inset map. Perhaps there is a better way to do this, but the way I've found is to use GMT's mapproject to calculate the height given the width. PyGMT doesn't expose mapproject so I have to use Session.call_module(), see the code below.

The problem is that doing this results in strange behaviour, dependent on where exactly the mapproject call is issued. I have found a way to arrange the code to get the figure I want, but it is very fragile, and the inset options do not work correctly.

Full code that generated the error

This is the "working" sample:

import pygmt
from pygmt.clib import Session
from pygmt.helpers import GMTTempFile


def mapproject_height(region, projection):
    with GMTTempFile() as tmp:
        args = f"-R{region} -J{projection} -Wh ->{tmp.name}"
        with Session() as lib:
            lib.call_module("mapproject", args)
        h = tmp.read().strip()

    return h


region = [-51, -48, -24, -23]
inset_region = [-80, -28, -43, 0]
inset_w = 2

fig = pygmt.Figure()

grid = pygmt.datasets.load_earth_relief(resolution="03s", region=region)
fig.grdimage(grid=grid, projection="M15c", frame=True, cmap="world")

fig.colorbar(frame=["x", "y+lm"], position="JBC+o0c/0.6c+w4.8c")

inset_h = mapproject_height("80W/28W/43S/0S", "M2c")
with fig.inset(position=f"jBL+w{inset_w}c/{inset_h}c+o0.1c", box="+pblack"):
    fig.coast(region=inset_region, projection=f"M{inset_w}c",
            land="gray", shorelines="thin", water="white")

fig.savefig("test.png")

Issues arise when trying to:

• Plot the colorbar after the inset---it is positioned relative to the inset rather than the main figure despite being outside the with statement.
• Move the inset to another corner (e.g. TL instead of BL). This results in the inset not being positioned correctly.
• Moving the mapproject_height call to the beginning, causes the outer figure to use the region boundaries from the inset for some reason.

Basically this only works if the inset is in the bottom-left, and plotted last, and the mapproject_height call is issued right before the with statement. This causes problems when I want to make more complex plots.

Full error message

n/a

System information

PyGMT information:
  version: v0.4.1
System information:
  python: 3.9.7 (default, Sep  3 2021, 12:37:55)  [Clang 12.0.5 (clang-1205.0.22.9)]
  executable: /Users/.../.local/share/virtualenvs/bed-mpQI0RL-/bin/python
  machine: macOS-11.6-x86_64-i386-64bit
Dependency information:
  numpy: 1.21.2
  pandas: 1.3.3
  xarray: 0.19.0
  netCDF4: 1.5.7
  packaging: 21.0
  ghostscript: 9.53.3
  gmt: 6.2.0
GMT library information:
  binary dir: /usr/local/Cellar/[email protected]/3.9.7/Frameworks/Python.framework/Versions/3.9/Resources/Python.app/Contents/MacOS
  cores: 8
  grid layout: rows
  library path: /usr/local/Cellar/gmt/6.2.0_1/lib/libgmt.dylib
  padding: 2
  plugin dir: /usr/local/Cellar/gmt/6.2.0_1/lib/gmt/plugins
  share dir: /usr/local/Cellar/gmt/6.2.0_1/share/gmt
  version: 6.2.0

[weiji14]

Hi @timrlaw, thanks for the detailed bug report! I haven't had time to fully reproduce this, but just to try and answer some of your questions quickly:

Essentially what I'm trying to do is set the dimensions of an inset figure to match the dimensions of the inset map.

Have you tried using fig.coast(..., projection="M?") to automatically determine the inset map's dimensions? See e.g. https://forum.generic-mapping-tools.org/t/hi-guys-would-you-help-me-find-that-why-my-inset-pic-is-not-complete-thanks/1675/2.

• Plot the colorbar after the inset---it is positioned relative to the inset rather than the main figure despite being outside the with statement.

• Move the inset to another corner (e.g. TL instead of BL). This results in the inset not being positioned correctly.

• Moving the mapproject_height call to the beginning, causes the outer figure to use the region boundaries from the inset for some reason.

Basically this only works if the inset is in the bottom-left, and plotted last, and the mapproject_height call is issued right before the with statement. This causes problems when I want to make more complex plots.

Having the colorbar is positioned relative to the inset map outside the with statement doesn't sound good. Could you upload some examples of how it looks? E.g. when setting the inset to be aligned at TL and BL.

[timrlaw]

Hi @weiji14

Sure thing, here are some pictures of how it looks. Starting from the code sample I put above, the image is "as desired":

If I move the colorbar call to be after the inset:

with fig.inset(...):
    ....
    
fig.colorbar(...)

If change "BL" on the inset to "TL", then I get this (certainly not in the top-left):

If I move mapproject_height(...) to the line right after inset_w = 2, I get this (on inspection, the outer plot is now covering the region specified for the inset plot, but only the zoomed in relief data is loaded, hence why it is so tiny):

I tried using "M?" on the inset fig.coast(...), but that didn't work either, I get the following:

I would also note that if I manually set inset_h to the literal value that mapproject_height gives, then everything works correctly. It is evidently caused by the Session.call_module().

[weiji14]

Ok, I figured it out. You need to put your mapproject_height call before doing fig = pygmt.Figure(), otherwise it messes up the GMT modern mode's hidden region and projection setting. I.e. do this:

region = [-51, -48, -24, -23]
inset_region = [-80, -28, -43, 0]
inset_w = 2
inset_h = mapproject_height(region="80W/28W/43S/0S", projection="M2c")  # put this before starting a new Figure!!

fig = pygmt.Figure()

grid = pygmt.datasets.load_earth_relief(resolution="03s", region=region)
fig.grdimage(grid=grid, projection="M15c", frame=True, cmap="world")
fig.colorbar(frame=["x", "y+lm"], position="JBC+o0c/0.6c+w4.8c")

with fig.inset(position=f"jTL+w{inset_w}c/{inset_h}c+o0.1c", box="+pblack"):
    fig.coast(
        region=inset_region,
        projection="M?",
        land="gray",
        shorelines="thin",
        water="white",
    )
fig.savefig("test.png")

produces

With this, you should be able to put the fig.colorbar line anywhere outside fig.inset and it should be positioned relative to the main figure.

This does seem similar to the global/local session problem @seisman mentioned in #733 (comment). What I don't get is that mapproject seems like a data processing function rather than a plotting function, so why should it affect the hidden GMT modern mode region and projection settings? Granted, I don't have a full understanding on how mapproject works, but it does seem like there is a bug here somewhere.

[maxrjones]

This does seem similar to the global/local session problem @seisman mentioned in #733 (comment). What I don't get is that mapproject seems like a data processing function rather than a plotting function, so why should it affect the hidden GMT modern mode region and projection settings? Granted, I don't have a full understanding on how mapproject works, but it does seem like there is a bug here somewhere.

There is no distinction between plotting and processing history for region and projection in GMT. I think this is a feature rather than a bug. For example, consider GMT's example 15. The region only needs to be defined in the first call (happens to be a processing module) and is then used throughout the remaining processing and plotting calls, unless specifically overwrote.

[timrlaw]

@weiji14, that code gives me the same result as my fourth picture above. I have tried with the version of PyGMT in PyPi and building from source. Are you perhaps using a development version of GMT? Could you please paste your pygmt.show_versions()?

@meghanrjones, that makes some sense. If I put a second mapproject_height call right after fig = pygmt.Figure() passing the zoomed in region coordinates (but discard the result of the function), then the main plot is the correct size. However the inset dimensions are now broken. I also note that the load_earth_relief call takes a region argument and this does not appear to reset the internal region.

This is quite confusing as a user, how do I know which PyGMT functions mutate the internal GMT state and which do not? Maybe I can make this work by scattering extra mapproject calls to reset the region before and after manipulating the inset, but this is contrary to how I would expect the with statement to work (i.e. the inset scope is not respected).

I have never used regular GMT, only PyGMT, so I have no mental model for how GMT works. However the API for PyGMT gives the impression that there is no internal state---based on how matplotlib and similar libraries work I would not expect calls made prior to fig = pygmt.Figure() to affect the plot. If this is indeed intended behaviour then perhaps it is worth clarifying somewhere prominent in the documentation.

Thanks both for your help.

[maxrjones]

Thanks for all your detailed reports. It is always helpful to get user feedback and non-GMT users provide an important perspective.

The improper setting of the inset dimensions may be an upstream bug; you shouldn't need to use mapproject to get this to work. I have opened an issue at GenericMappingTools/gmt#5875 and will post back with the relevant results for PyGMT.

Regarding the state changes, I agree that calls before creating a figure instance should not affect the plot. I found that @weiji14's example worked well in the development branch but failed in v0.4.1, as you noticed. Edit: Actually, I am getting the same as the fourth figure with both the dev branch and v0.4.1. From GMT's perspective, all settings from a session are accessible to the figure level. Perhaps Wei Ji or others have information about whether this is intended to be the case for PyGMT as well.

[maxrjones]

Independent of the issues you're having with the region setting, using mapproject is the correct way to get the inset size to agree with the region and projection. We'll look into simplifying this upstream.

[timrlaw]

Thanks @meghanrjones. For reference, my use case is automatically generating plots to be displayed on a website. Hence why I cannot manually fiddle the inset height---it needs to be fully programmable.

[weiji14]

Really sorry for the confusion @timrlaw, I posted the wrong code 😅, forgot to put a region argument in the fig.grdimage call (which means it was picking up the region from mapproject). Try this instead:

region = [-51, -48, -24, -23]
inset_region = [-80, -28, -43, 0]
inset_w = 2
inset_h = mapproject_height(region="80W/28W/43S/0S", projection="M2c")  # put this before starting a new Figure!!

fig = pygmt.Figure()

grid = pygmt.datasets.load_earth_relief(resolution="03s", region=region)
fig.grdimage(grid=grid, region=region, projection="M15c", frame=True, cmap="world")  # need to reset the region here
fig.colorbar(frame=["x", "y+lm"], position="JBC+o0c/0.6c+w4.8c")

with fig.inset(position=f"jTL+w{inset_w}c/{inset_h}c+o0.1c", box="+pblack"):
    fig.coast(
        region=inset_region,
        projection="M?",
        land="gray",
        shorelines="thin",
        water="white",
    )
fig.savefig("test.png")

For reference, this is my pygmt.show_versions(), using try-gmt

PyGMT information:
  version: v0.4.1
System information:
  python: 3.9.7 | packaged by conda-forge | (default, Sep 23 2021, 07:28:37)  [GCC 9.4.0]
  executable: /srv/conda/envs/notebook/bin/python
  machine: Linux-5.4.129+-x86_64-with-glibc2.27
Dependency information:
  numpy: 1.21.2
  pandas: 1.3.3
  xarray: 0.19.0
  netCDF4: 1.5.7
  packaging: 20.4
  ghostscript: 9.54.0
  gmt: 6.2.0
GMT library information:
  binary dir: /srv/conda/envs/notebook/bin
  cores: 6
  grid layout: rows
  library path: /srv/conda/envs/notebook/lib/libgmt.so
  padding: 2
  plugin dir: /srv/conda/envs/notebook/lib/gmt/plugins
  share dir: /srv/conda/envs/notebook/share/gmt
  version: 6.2.0

This is quite confusing as a user, how do I know which PyGMT functions mutate the internal GMT state and which do not? Maybe I can make this work by scattering extra mapproject calls to reset the region before and after manipulating the inset, but this is contrary to how I would expect the with statement to work (i.e. the inset scope is not respected).

I have never used regular GMT, only PyGMT, so I have no mental model for how GMT works. However the API for PyGMT gives the impression that there is no internal state---based on how matplotlib and similar libraries work I would not expect calls made prior to fig = pygmt.Figure() to affect the plot. If this is indeed intended behaviour then perhaps it is worth clarifying somewhere prominent in the documentation.

Yes this is good feedback, and is essentially what we have been discussing at #1282 on having unintended behaviour arising from using a single GMT session for everything in PyGMT. It is a tricky problem which Leo appears to have debated on early in PyGMT's history (see #327 (comment)).

@meghanrjones, do you think it would be better to let fig = pygmt.Figure() reset the internal state (i.e. clear what is in the .gmt/sessions/gmt_session.12345 folder)? Perhaps through an option like pygmt.Figure(reset=True)? At the very least, the region and projection could be reset somehow so people don't get maps with the wrong geographic extent plotted. Or we could add a note somewhere in the class Figure docstring to mention this caveat that there could be stuff affecting the Figure's state somewhere here:

pygmt/pygmt/figure.py

Lines 47 to 56 in 65995f6

	class Figure:
	"""
	A GMT figure to handle all plotting.
	Use the plotting methods of this class to add elements to the figure. You
	can preview the figure using :meth:`pygmt.Figure.show` and save the figure
	to a file using :meth:`pygmt.Figure.savefig`.
	Unlike traditional GMT figures, no figure file is generated until you call
	:meth:`pygmt.Figure.savefig` or :meth:`pygmt.Figure.psconvert`.

[timrlaw]

Thanks @weiji14, that works a treat, I think I can do everything I need now. As long as I do all my mapprojecting up front before fig = pygmt.Figure(), and then make sure to set the region on the first drawing call, changing the order of the subsequent elements of the plot all seems to work fine.

I will let you decide whether or not to close this issue now or keep it around pending possible changes.

Thanks again both for your help.

[maxrjones]

@meghanrjones, do you think it would be better to let fig = pygmt.Figure() reset the internal state (i.e. clear what is in the .gmt/sessions/gmt_session.12345 folder)? Perhaps through an option like pygmt.Figure(reset=True)? At the very least, the region and projection could be reset somehow so people don't get maps with the wrong geographic extent plotted. Or we could add a note somewhere in the class Figure docstring to mention this caveat that there could be stuff affecting the Figure's state somewhere here:

I like your idea of providing an option for pygmt.Figure(reset=True). We could still have the default as False to maintain backwards compatibility. We could either end the session and begin a new session if reset=True or have an upstream feature request for gmt clear history and then call that if reset=True. I would lean towards the second option. What do you think?

In the meantime, this seems to be a recurring tricky point with the history of parameters when dealing with figures, subplots, and insets. I think we should add some documentation to clarify these points and will open a separate issue to track that.