-
Notifications
You must be signed in to change notification settings - Fork 225
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
grdimage does not work correctly when the redundant 360 E latitude is not included #3331
Comments
@MarkWieczorek thanks for reporting this problem! I tested your code, and there seems to be an OS-dependency:
import numpy as np
import xarray as xr
import pygmt
# grid that does not include redundant 360 E
grid2 = np.random.rand(181, 360)
lat2 = np.arange(90, -91, -1)
lon2 = np.arange(0, 360, 1)
da2 = xr.DataArray(grid2, coords=[("lat", lat2), ("lon", lon2)])
da2.gmt.registration = 0
da2.gmt.gtype = 1
figure = pygmt.Figure()
figure.grdimage(da2, projection="H10c", region="g")
figure.coast(shorelines="1/1p,blue")
figure.shift_origin(xshift="+w+1c")
figure.grdimage(da2, projection="H10c", region="d")
figure.coast(shorelines="1/1p,blue")
figure.shift_origin(xshift="-w-1c", yshift="-h-1c")
figure.grdimage(da2, projection="H60/10c", region="g")
figure.coast(shorelines="1/1p,blue")
figure.shift_origin(xshift="+w+1c")
figure.grdimage(da2, projection="H60/10c", region="d")
figure.coast(shorelines="1/1p,blue")
figure.show()
|
Same problem on my intel macOS machine. |
I can reproduce your issue. If I change I think it's an upstream bug. Will see if I can fix it later. |
Additionally, I found there is a problem for the case the longitude is given in the range 0°-360° E, but not for -180°-180° E. import xarray as xr
import pygmt
grid = np.random.rand(180, 360)
lat = np.arange(90, -90, -1)
lon2 = np.arange(0, 360, 1)
lon3 = np.arange(-180, 180, 1)
da2 = xr.DataArray(grid, coords=[("lat", lat), ("lon", lon2)])
da2.gmt.registration = 0
da2.gmt.gtype = 1
da3 = xr.DataArray(grid, coords=[("lat", lat), ("lon", lon3)])
da3.gmt.registration = 0
da3.gmt.gtype = 1
for da in [da2, da3]:
figure = pygmt.Figure()
figure.grdimage(da, projection="N10c", region="g")
figure.coast(shorelines="1/1p,blue", frame=["af", "+tN10c | g"])
figure.shift_origin(xshift="+w+1c")
figure.grdimage(da, projection="N10c", region="d")
figure.coast(shorelines="1/1p,blue", frame=["af", "+tN10c | d"])
figure.shift_origin(xshift="-w-1c", yshift="-h-2c")
figure.grdimage(da, projection="N0/10c", region="g")
figure.coast(shorelines="1/1p,blue", frame=["af", "+tN0/10c | g"])
figure.shift_origin(xshift="+w+1c")
figure.grdimage(da, projection="N180/10c", region="d")
figure.coast(shorelines="1/1p,blue", frame=["af", "+tN180/10c | d"])
figure.show()
|
The error is caused by an upstream bug, which is fixed in GenericMappingTools/gmt#8554. |
The upstream PR GenericMappingTools/gmt#8554 has been merged. With the GMT dev version installed, the following test shows now it works well for any centering longtiude:
A similar test is added in #3358. Please note that, if you run the example in #3331 (comment), save the two images into two PNG files and then compare them closely, you will still notice minor differences at the map edges and the 180° longitude line. These differences are mostly caused by paddings and they will more likely be fixed in #3099. |
Here is a more comprehensive test for 4 different grids:
These 4 grids are plotted with different projection central longitude. import numpy as np
import xarray as xr
import pygmt
from pygmt.datasets import load_earth_relief
# Global grid in the longitude [-180, 180] range with redundant longitude at 180/-180
da1 = load_earth_relief(region=[-180, 180, -90, 90])
assert da1.gmt.registration == 0
assert da1.gmt.gtype == 1
assert da1.shape == (181, 361)
assert da1.lon.values.min() == -180.0
assert da1.lon.values.max() == 180.0
# Global grid in the longitude [0, 360] range with redundant longitude at 0/360
da2 = load_earth_relief(region=[0, 360, -90, 90])
assert da2.gmt.registration == 0
assert da2.gmt.gtype == 1
assert da2.shape == (181, 361)
assert da2.lon.values.min() == 0.0
assert da2.lon.values.max() == 360.0
# Global grid in the longitude [-180, 180] range without redundant longitude at 180/-180
da3 = da1[:, 0:360]
da3.gmt.registration, da3.gmt.gtype = 0, 1
assert da3.shape == (181, 360)
assert da3.lon.values.min() == -180.0
assert da3.lon.values.max() == 179.0
# Global grid in the longitude [0, 360] range without redundant longitude at 0/360
da4 = da2[:, 0:360]
da4.gmt.registration, da4.gmt.gtype = 0, 1
assert da4.shape == (181, 360)
assert da4.lon.values.min() == 0.0
assert da4.lon.values.max() == 359.0
fig = pygmt.Figure()
for lon in np.arange(5, 360, 30):
kwdict = {
"projection": f"W{lon}/10c",
"region": "g",
"frame": f"+t{lon:.2f}"
}
fig.grdimage(da1, **kwdict)
fig.shift_origin(xshift="w+2c")
fig.grdimage(da2, **kwdict)
fig.shift_origin(xshift="w+2c")
fig.grdimage(da3, **kwdict)
fig.shift_origin(xshift="w+2c")
fig.grdimage(da4, **kwdict)
fig.shift_origin(xshift="-3w-6c", yshift="h+2c")
fig.show() Each row contains results for the 4 images with a specific central longitude (shown as title). As can be seen, the bug exists under the following conditions: (1) grid longitude is in the [0, 360] range; (2) no redundant longitude at 0/360; (3) central longitude is < 180. |
Description of the problem
grdimage works fine for global grids that include the redundant 0 and 360 E longitudes. However, if you create a grid that doesn't include the redundant 360 E, the returned image is all black.
I don't know if this worked before or not. For my purposes, it would be useful to use grdimage with global arrays that exclude 360 E and perhaps 90 S. The reason for this is that the arrays used for creating spherical harmonic transforms often don't require these values.
Minimal Complete Verifiable Example
Full error message
System information
The text was updated successfully, but these errors were encountered: