Some geographic mapping issues

[chrishavlin]

Starting this issue to document some issues with yt's geographic plotting functionality.

Test data and environment

Links to the data used below, both are seismic models (internal_geographic in yt jargon) the first is a global model (longitude from -180 to 180, latitude -90 to 90), the second covers only the western US.

• https://ds.iris.edu/files/products/emc/emc-files/GYPSUM_percent.nc
• https://ds.iris.edu/files/products/emc/emc-files/wUS-SH-2010_percent.nc

And here's a function for loading that data in as a internal_geographic uniform grid datasets: https://gist.github.com/chrishavlin/aaf887a5f36eb5f6a7f7ae45f936c2fa

The following used a conda environment with python 3.9.10, cartopy 0.20.3 , geos 3.11.0 , proj 9.0.1 yt 4.1.dev0 on unbuntu 18.04.

Problem 1: extents

I know there's been some issues here previously (#3705, #3708), but the current state does not do so well...

The default behavior for a global dataset is good:

ds = get_internal_IRIS("internal_geographic_global")
c = ds.domain_center
c[0] = 100
p = yt.SlicePlot(ds,'depth','dvs', center=c)
p.set_mpl_projection('Robinson')
p.set_log('dvs',False)
p.set_zlim('dvs',-8, 8)
p.set_cmap('dvs','magma_r')

p._setup_plots()
carto_ax = p.plots['dvs'].axes
carto_ax.add_feature(cfeature.NaturalEarthFeature(
        'cultural', 'admin_0_countries', '10m',
        edgecolor='black', facecolor='none'))
p.show()

But any subsets are VERY wrong:

c = ds.domain_center
c[0] = 100
p = yt.SlicePlot(ds,'depth','dvs', center=c, width=[100, 60], origin="native")
p.set_log('dvs',False)
p.set_zlim('dvs',-8, 8)
p.set_cmap('dvs','magma_r')
p._setup_plots()
carto_ax = p.plots['dvs'].axes
carto_ax.add_feature(cfeature.NaturalEarthFeature(
        'cultural', 'admin_0_countries', '10m',
        edgecolor='black', facecolor='none'))
p.show()

This happens because the extent gets set to None for cartopy GeoAxes (happens here) and so the image slice gets stretched to fill the globe. Commenting out that extent = None line results in a correct image:

Note that the GeoAxes extent is still the whole globe and can be zoomed in with a call the axis' set_extent method:

carto_ax.set_extent((-50, 50, -30, 30))
p.show()

For a dataset that does not cover the whole globe, it immediately gets stretched:

ds = get_internal_IRIS("internal_geographic_partial")

c = ds.domain_center
c[0] = 100
p = yt.SlicePlot(ds,'depth','dvs', center=c)
p.set_mpl_projection('Robinson')
p.set_log('dvs',False)
p.set_zlim('dvs',-8, 8)
p.set_cmap('dvs','magma_r')

p._setup_plots()
carto_ax = p.plots['dvs'].axes
carto_ax.add_feature(cfeature.NaturalEarthFeature(
        'cultural', 'admin_0_countries', '10m',
        edgecolor='black', facecolor='none'))
p.show()

Again commenting out that extent=None line, it gets correctly placed in the globe:

Given the cartopy incompatibility issues (#3705), simply removing the extents = None is obviously not an option. It would be possible only set extents=None when the extent covers the globe, but I'll have to check if that will expose previous issues...

Problem 2: rotating a view

This is likely related, but worth mentioning separately: rotating the view is not straightforward. For reference, here's the default plot for the global data again:

Trying to rotate the view with yt's center argument does rotate the selection correctly but since this doesn't register anything with the GeoAxes, subsequent plots do not respect the rotation:

c = ds.domain_center
c[0] = 100
c[2] = 50
p = yt.SlicePlot(ds,'depth','dvs', center=c)
p.set_log('dvs',False)
p.set_zlim('dvs',-8, 8)
p.set_cmap('dvs','magma_r')

p._setup_plots()
carto_ax = p.plots['dvs'].axes
carto_ax.add_feature(cfeature.NaturalEarthFeature(
        'cultural', 'admin_0_countries', '10m',
        edgecolor='black', facecolor='none'))
p.show()

If you're aware of this behavior, you can avoid it by adjusting the projection instead of the center argument:

c = ds.domain_center
c[0] = 100
p = yt.SlicePlot(ds,'depth','dvs', center=c)
p.set_mpl_projection(('Mollweide', (), {'central_longitude': 50}))
p.set_log('dvs',False)
p.set_zlim('dvs',-8, 8)
p.set_cmap('dvs','magma_r')

p._setup_plots()
carto_ax = p.plots['dvs'].axes
carto_ax.add_feature(cfeature.NaturalEarthFeature(
        'cultural', 'admin_0_countries', '10m',
        edgecolor='black', facecolor='none'))
p.show()

But it'd be nice if the default projection takes into account the SlicePlot center argument.

[neutrinoceros]

Given the cartopy incompatibility issues (#3705), simply removing the extents = None is obviously not an option.

I wouldn't be so sure actually. We have very limited venue to test integration with cartopy (at the moment we only test it on windows + conda) because it's so hard to setup, and I don't think we do any kind of automated image testing for it, which might explain how I got it so bad when I wrote this line. Maybe this line could be conditional (depending on cartopy's version), or maybe we could just clean it up and only claim support with the latest version of cartopy (which is the only one we test for at any point anyway). Anyway my point is: don't treat this line as sacred. I may only have been smashing my keyboard until CI got green again.

And thank you very much for all the work you've done to report these problem !

[matthewturk]

I'm not totally sure that I understand -- if I am reading it right, the problem with having it not call set_extent at all was for subsets of the globe, right? But now that seems to work?

And if I am understanding, I think that if we have a minimum cartopy version set, we should just ... be OK? So maybe just avoiding calling set_extent is the solution? Here, I find that removing the line works for both full and subset with no hangs.

[chrishavlin]

Update from some discussion between me and Matt -- we decided that we probably do want to call set_extent.

I just pushed up a draft ( #4040 ) that only calls it if the plot is not a global plot (within the tolerance of the transform). I think that might do a good job... going to still look into potentially setting a central_longitude for the projection when possible (and if feasible).