Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Ensure plotting xarray grids with different central meridians work on some projections #560

Merged
merged 19 commits into from
Sep 11, 2020

Conversation

weiji14
Copy link
Member

@weiji14 weiji14 commented Aug 10, 2020

Description of proposed changes

Adds tests to make sure that plotting xarray.DataArray grids at different central meridians work properly for different projection types. See #515 (comment) for background context. Note that this depends on GMT >= 6.1.1 and won't work yet.

Reverts the runfirst workaround required in #476, reverts #531 grdtrack gallery example crash on macOS due to a problematic meridian setting.

Note that this PR may depend on #555 ✔️

Addresses #390, Fixes #515

Reminders

  • Run make format and make check to make sure the code follows the style guide.
  • Add tests for new features or tests that would have caught the bug that you're fixing.
  • Add new public functions/methods/classes to doc/api/index.rst.
  • Write detailed docstrings for all functions/methods.
  • If adding new functionality, add an example to docstrings or tutorials.

@weiji14 weiji14 added the bug Something isn't working label Aug 10, 2020
@weiji14 weiji14 added this to the 0.2.x milestone Sep 2, 2020
@weiji14
Copy link
Member Author

weiji14 commented Sep 2, 2020

As a quick test, I've reverted #531, and the grdtrack gallery example works fine again! See https://pygmt-git-test-centralmeridians.gmt.now.sh/gallery/grid/track_sampling.html, and compare it with screenshots at #515 (comment).

Asia Pacific centred grdtrack example

@weiji14 weiji14 changed the title WIP Ensure plotting xarray grids with different central meridians work Ensure plotting xarray grids with different central meridians work Sep 3, 2020
@weiji14
Copy link
Member Author

weiji14 commented Sep 3, 2020

Ok, so plotting at different central meridians doesn't crash anymore (and I'm quite happy with that)! However, I'm not sure how I feel about the bright green pixels on the rightmost edge. The below image is for gridline registered, geographic type xarray.DataArray grids. Top row are NetCDF grids, middle and bottom rows are from xarray.DataArray grids:

test

Code from #390 (comment), adapted slightly.

import numpy as np
import pygmt
import xarray as xr

# Creata a data array in gridline coordinates of sin(lon) * cos(lat)
interval = 10
lat = np.arange(90, -90 - interval, -interval)
lon = np.arange(0, 360 + interval, interval)
longrid, latgrid = np.meshgrid(lon, lat)
data = np.sin(np.deg2rad(longrid)) * np.cos(np.deg2rad(latgrid))

# create a DataArray
dataarray = xr.DataArray(
    data,
    coords=[
        ("latitude", lat, {"units": "degrees_north"}),
        ("longitude", lon, {"units": "degrees_east"}),
    ],
    attrs={"actual_range": [-1, 1]},
)
assert dataarray.gmt.registration == 0  # Gridline registration
dataarray.gmt.gtype = 1  # Geographic
dataarray.to_dataset(name="dataarray").to_netcdf("test.grd")
print(pygmt.grdinfo("test.grd"))

# create projected images using a cylindrical and mollweide projection.
fig = pygmt.Figure()

regframeargs = {"region": "g", "frame": "a90f30g30"}
# Lower Left
fig.grdimage(dataarray, projection="Q0/0/6i", **regframeargs)
# Lower Right
fig.grdimage(dataarray, projection="Q180/0/6i", X="7i", **regframeargs)
# Middle Left
fig.grdimage(dataarray, projection="W0/6i", X="-7i", Y="4i", **regframeargs)
# Middle Right
fig.grdimage(dataarray, projection="W180/6i", X="7i", **regframeargs)
# Top Left
fig.grdimage("test.grd", projection="Q0/0/6i", X="-7i", Y="4i", **regframeargs)
# Top Right
fig.grdimage("test.grd", projection="Q180/0/6i", X="7i", **regframeargs)

fig.savefig("test.png")
fig.show()

Maybe we just need to add another column of pixels on the Eastern side? Will look into it after lunch.

Testing meridians 0, 33, 120 and 180;
projections H, Q and W. This requires
modifications to the check_figures_equal
decorator function so that 
pytest.mark.parametrize can work.
@weiji14
Copy link
Member Author

weiji14 commented Sep 4, 2020

The Cylindrical Equidistant ('Q') style projection tests fail for meridians centred on 33, 120 and 180 (works for 0 deg). Below is for the 120deg test, the difference seems to be more of a pixelation than anything. Quite strange.

Using netcdf Using xarray diff
test_grdimage_different_central_meridians_and_projections Q-120-png -expected test_grdimage_different_central_meridians_and_projections Q-120-png test_grdimage_different_central_meridians_and_projections Q-120-png -failed-diff

Other tested projections Hammer (H) and Mollweide (W) work fine for 0, 33, 120 and 180 deg centred maps.

Split up tests into projections that only take a
central meridian argument (e.g. H, W); and
those that also take a standard parallel
argument (e.g. Q and T). Reduced test scope
to only test on longitudes 0, 123 and 180.
Copy link
Member Author

@weiji14 weiji14 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For anyone following (cc @MarkWieczorek), this PR is basically the main one we'll like to merge in for getting a PyGMT v0.2.0 out. Without getting too technical, the recent GMT 6.1.1 release has solved a lot of things for us, such as fatal crashes and black plots when using xarray DataArray grids. However, the tests in this PR show that a few gotchas remain, particularly with the Cylindrical Equidistant (Q) projection issue mentioned at #390.

How important is it that we get this perfected? I'm entertaining the idea of getting a PyGMT v0.2.0 release out as soon as possible instead of waiting for a GMT 6.1.2 release (or a GMT 6.2.0 release in November). I'm not sure if there's much we can do on the PyGMT side (besides testing the GMT C API), and we can always do a PyGMT v0.2.1 release once these problems are resolved in upstream GMT (by mid to late September?).

Comment on lines 96 to 106
@pytest.mark.parametrize("lon0", [0, 123, 180])
@pytest.mark.parametrize("proj_type", ["H", "W"])
def test_grdimage_central_meridians(grid, proj_type, lon0, fig_ref, fig_test):
"""
Test that plotting a grid centred at different longitudes/meridians work.
Test that plotting a grid with different central meridians (lon0) using
Hammer (H) and Mollweide (W) projection systems work.
"""
fig_ref.grdimage("@earth_relief_01d_g", projection="W120/15c", cmap="geo")
fig_test.grdimage(grid, projection="W120/15c", cmap="geo")
fig_ref.grdimage(
"@earth_relief_01d_g", projection=f"{proj_type}{lon0}/15c", cmap="geo"
)
fig_test.grdimage(grid, projection=f"{proj_type}{lon0}/15c", cmap="geo")
Copy link
Member Author

@weiji14 weiji14 Sep 6, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This test checks that #515 is fixed (and it is!), specifically, that PyGMT can plot xarray grids using the Hammer (H) (and Mollweide (W)) projection system on central longitudes 0, 123 and 180deg.

Comment on lines 109 to 124
@check_figures_equal()
@pytest.mark.parametrize("lat0", [0, 30])
@pytest.mark.parametrize("lon0", [0, 123, 180])
@pytest.mark.parametrize("proj_type", ["Q", "S"])
def test_grdimage_central_meridians_and_standard_parallels(
grid, proj_type, lon0, lat0, fig_ref, fig_test
):
"""
Test that plotting a grid with different central meridians (lon0) and
standard_parallels (lat0) using Cylindrical Equidistant (Q) and General
Stereographic (S) projection systems work.
"""
fig_ref.grdimage(
"@earth_relief_01d_g", projection=f"{proj_type}{lon0}/{lat0}/15c", cmap="geo"
)
fig_test.grdimage(grid, projection=f"{proj_type}{lon0}/{lat0}/15c", cmap="geo")
Copy link
Member Author

@weiji14 weiji14 Sep 6, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This check is meant to resolve @MarkWieczorek's longstanding issue at #390 (we're getting there). Specifically, we check plotting xarray grids on different central meridians (lon0) at 0, 123 and 180deg, and different standard parallels (lat0) 0, 30deg for Cylindrical Equidistant (Q) and General Stereographic (S) projections systems. These plots are compared to a baseline/reference plot using a NetCDF grid directly.

Good news: With GMT 6.1.1, plotting these grids don't crash or give a black image anymore. Also for a few lon0/lat0 combinations e.g. (Q0/0, Q0/30, S123/0, S180/0), the plots are perfect.

Bad news: Certain lon0/lat0 combinations plot in xarray, but aren't exactly the same as the NetCDF plots. There's an RMS difference of ~25 for the Cylindrical Equidistant plots Q123/0, Q123/30, Q180/0, Q180/30 (see my comment at #560 (comment) for a visual diff). The General Stereographic plots also show minor RMS differences <2 for S0/0, S0/30, S123/30, S180/30, with the diff appearing as a strip along 0deg longitude (i.e. Greenwich, see #560 (comment) for a visual diff).

I've also tested combinations of GMT_IN|GMT_IS_REFERENCE (the current setting), GMT_IN (xref #517), and GMT_IN|GMT_IS_DUPLICATE, but the results are the same. Paul's comment at GenericMappingTools/gmt#3829 (comment) does suggest that there might be outstanding issues to resolve here.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for taking to the time to look into this. My attitude is that if you think you have done all you can with GMT 6.1.1, then you should release pygmt 0.2 as soon as possible, even if there are still some bugs that might be upstream.

pygmt 0.1.2 is completely unusable for me (for plotting global projected data with xarray grids), and as far as I can tell with my tests from pygmt/master, things appear to be ok (i.e., no kernel crashes, and the images appear to be fine without detailed inspection). Nevertheless, the problem with the last column on the right of the plot being incorrect is a clear bug and will need to resolved in the following release.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It seems there is nothing we can do in PyGMT to fix these bugs. We could do a 0.2.0 release soon.

@weiji14 weiji14 mentioned this pull request Sep 6, 2020
11 tasks
@vercel vercel bot temporarily deployed to Preview September 8, 2020 03:04 Inactive
Also set RMS tolerance to 1.5 to handle small different along Greenwich Meridian for General Stereographic (S) plots.
Comment on lines +112 to +120
# Cylindrical Equidistant (Q) projections plotted with xarray and NetCDF grids
# are still slightly different with an RMS error of 25, see issue at
# https://github.com/GenericMappingTools/pygmt/issues/390
# TO-DO remove tol=1.5 and pytest.mark.xfail once bug is solved in upstream GMT
@check_figures_equal(tol=1.5)
@pytest.mark.parametrize("lat0", [0, 30])
@pytest.mark.parametrize("lon0", [0, 123, 180])
@pytest.mark.parametrize("proj_type", [pytest.param("Q", marks=pytest.mark.xfail), "S"])
def test_grdimage_central_meridians_and_standard_parallels(grid, proj_type, lon0, lat0):
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Using xfail (again) to temporarily allow the Cylindrical Equidistant (Q) tests to pass (instructions at https://docs.pytest.org/en/latest/skipping.html#skip-xfail-with-parametrize). Note that there are **2 xpass**es (because not all projections are wrong), and 4 xfails.

Also, I've set the RMS tolerance to 1.5 to allow the General Stereographic (S) projection tests to pass, for reference, here are the plot diffs (along the Greenwich Meridian):

NetCDF (expected) Xarray (test) Diff
test_grdimage_central_meridians_and_standard_parallels S-0-0-png -expected test_grdimage_central_meridians_and_standard_parallels S-0-0-png test_grdimage_central_meridians_and_standard_parallels S-0-0-png -failed-diff

@vercel vercel bot temporarily deployed to Preview September 11, 2020 03:11 Inactive
@weiji14 weiji14 force-pushed the test/central_meridians branch from cbaed26 to fea219f Compare September 11, 2020 03:13
@vercel vercel bot temporarily deployed to Preview September 11, 2020 03:13 Inactive
@weiji14 weiji14 changed the title Ensure plotting xarray grids with different central meridians work Ensure plotting xarray grids with different central meridians work on some projections Sep 11, 2020
@weiji14 weiji14 force-pushed the test/central_meridians branch from fea219f to b3e063a Compare September 11, 2020 03:17
@weiji14 weiji14 marked this pull request as ready for review September 11, 2020 03:27
Copy link
Member

@seisman seisman left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It seems you're doing two things in one PR. Can it be two separate PRs?

@weiji14
Copy link
Member Author

weiji14 commented Sep 11, 2020

It seems you're doing two things in one PR. Can it be two separate PRs?

Yep, done at #600, please review that one before we merge in this one.

Copy link
Member

@seisman seisman left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks good, just a few minor comments.

pygmt/tests/test_grdimage.py Outdated Show resolved Hide resolved
pygmt/tests/test_grdimage.py Outdated Show resolved Hide resolved
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug Something isn't working
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Global earth relief data loaded as xarray doesn't work if projection center is not 0
3 participants