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

Matrix as grid with changing central meridian doesn't work well for gridline grids #3844

Closed
seisman opened this issue Aug 5, 2020 · 15 comments
Labels
bug Something isn't working

Comments

@seisman
Copy link
Member

seisman commented Aug 5, 2020

Full script that generated the error

The script generates 4 maps, for pixel or gridline grids passing as grid or matrix:

import pygmt
from pygmt.datasets import load_earth_relief
import xarray as xr

# pixel grid
fig = pygmt.Figure()
fig.grdimage("@earth_relief_01d_p", projection="W120/15c", cmap="geo")
fig.savefig("grid_p.png")

# pixel matrix
grid = load_earth_relief(registration="pixel")
fig = pygmt.Figure()
fig.grdimage(grid, projection="W120/15c", cmap="geo")
fig.savefig("matrix_p.png")


# gridline grid
fig = pygmt.Figure()
fig.grdimage("@earth_relief_01d_g", projection="W120/15c", cmap="geo")
fig.savefig("grid_g.png")

# gridline matrix
grid = load_earth_relief(registration="gridline")
fig = pygmt.Figure()
fig.grdimage(grid, projection="W120/15c", cmap="geo")
fig.savefig("matrix_g.png")

Use the command below to compare the two maps of pixel-registered grid/matrix. The RMS is zero. So passing grid or matrix give the same results:

gm compare -density 200 -highlight-color magenta -highlight-style assign -metric rmse -file diff_p.png matrix_p.png grid_p.png

Use the command below to compare the two maps of pixel-registered grid/matrix.

gm compare -density 200 -highlight-color magenta -highlight-style assign -metric rmse -file diff_g.png matrix_g.png grid_g.png

The RMS is non-zero. Here is the diff image. There are tiny difference along the x=180 line. Perhaps a 0.5 grid offset?

image

@seisman seisman added the bug Something isn't working label Aug 5, 2020
@PaulWessel
Copy link
Member

Will look tomorrow. I know the issue: When we rotate the central longitude the repeated columns at -180 and 180 need to collapse to one.

@PaulWessel
Copy link
Member

So 1-2 days go by, an automatic OS update happened over night, so GMT_LIBRARY_PATH was unset. However, I am missing another step before I can sy import pygmt, something about switching from base to pygmt or a setup step?

@seisman
Copy link
Member Author

seisman commented Aug 5, 2020

Yes, you need to activate your conda environment by conda activate pygmt.

@PaulWessel
Copy link
Member

??? My outputs look like this (with master):

matrix_g

@seisman
Copy link
Member Author

seisman commented Aug 5, 2020

The script generates 4 figures. Which one?

@PaulWessel
Copy link
Member

They are all like that (just a sliver near the boundary), I think this one was matrix_g.png

@seisman
Copy link
Member Author

seisman commented Aug 5, 2020

This makes no sense. Perhaps try pygmt.show_versions() and see if you're using the correct GMT or PyGMT version?

@PaulWessel
Copy link
Member

Makes perfect sense. There was an old gmt.history in the directory I started python:

GMT 6 Session common arguments shelf

BEGIN GMT 6.2.0
B 0g10
J M
JM M15c
JQ Q15/6i
R -115/-100/55/60
@l 1
END

and the python script latched onto that. SHould not pygmt start a new session and no gmt.history is used?

@seisman
Copy link
Member Author

seisman commented Aug 5, 2020

You mean the gmt begin -C?

Start this session with a clean slate: Any gmt.conf files in the usual search path directories are ignored [Default starts session with the prevailing user settings].

Does it also ignore gmt.history?

@PaulWessel
Copy link
Member

So back to the original issue: Since we are passing reference we simply use the grid. It gets projected, in gmt_grd_project there will be loops over input and output nodes (you remember this from last week). So there will be two sets of meridians taht are used (lon = -180 and lon = 180) which are exactly the same. So those should technically be given a 0.5 weight so that in averaging many nodes they should not count more. Normally (i.e., grdfiles) we do not see this (or possibly never noticed) because we read int the grid with rotation, so the repeated columns are at the two separate sides. Now they show up in the middle.

I will see if ti is possibly to insert the weight of 0.5.

@PaulWessel
Copy link
Member

If you start with begin -C then clearly we leak that through. Give it a test bu using my history file above and confirm the problem.

@PaulWessel
Copy link
Member

Testing a fix for this. Now it passes:
Normalized Absolute
============ ==========
Red: 0.0054440957 1.4
Green: 0.0035965496 0.9
Blue: 0.0052102560 1.3
Total: 0.0048207923 1.2

versus before it was

       Normalized    Absolute
      ============  ==========
 Red: 0.0076255501        1.9

Green: 0.0050802491 1.3
Blue: 0.0066723388 1.7
Total: 0.0065441590 1.7

All tests still pass so I am making a PR.

@PaulWessel
Copy link
Member

Since the history file is read in GMT_Create_Session, it has already been ingested when begin -C is parsed. Thus, I need to reset the history. I could not recreate the problem on the command line, so suspect it is interfering with the matri/grdimage stuff.

@PaulWessel
Copy link
Member

There should be no instances where an old gmt.history file should be used by a new modern script, -C or not. I think I will always reset history when gmt begin is called. It is different for the settings.

@seisman
Copy link
Member Author

seisman commented Aug 6, 2020

Closed by #3849.

@seisman seisman closed this as completed Aug 6, 2020
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

No branches or pull requests

2 participants