grdimage does not work correctly when the redundant 360 E latitude is not included

MarkWieczorek commented 3 months ago

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

import numpy as np
import xarray as xr
import pygmt

# grid that includes redundant 360 E
grid = np.random.rand(181, 361)
lat = np.arange(90, -91, -1)
lon = np.arange(0, 361, 1)

da = xr.DataArray(grid, coords=[('lat', lat), ('lon', lon)])
da.gmt.registration = 0
da.gmt.gtype = 1

figure = pygmt.Figure()
figure.grdimage(da, projection='W0/10c', region='g')
figure.show()

# 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

figure2 = pygmt.Figure()
figure2.grdimage(da2, projection='W0/10c', region='g')
figure2.show()

Full error message

There is no error. The image is created, but it is all black.

System information

PyGMT information:
  version: v0.12.0
System information:
  python: 3.12.3 | packaged by conda-forge | (main, Apr 15 2024, 18:38:13) [GCC 12.3.0]
  executable: /home/mrak/miniforge3/envs/py312/bin/python
  machine: Linux-6.9.7-200.fc40.x86_64-x86_64-with-glibc2.39
Dependency information:
  numpy: 1.26.4
  pandas: 2.2.2
  xarray: 2024.3.0
  netCDF4: 1.6.5
  packaging: 24.0
  contextily: None
  geopandas: None
  ipython: None
  rioxarray: None
  ghostscript: 10.03.0
GMT library information:
  binary version: 6.5.0
  cores: 20
  grid layout: rows
  image layout: 
  library path: /home/mrak/miniforge3/envs/py312/lib/libgmt.so
  padding: 2
  plugin dir: /home/mrak/miniforge3/envs/py312/lib/gmt/plugins
  share dir: /home/mrak/miniforge3/envs/py312/share/gmt
  version: 6.5.0

yvonnefroehlich commented 3 months ago

@MarkWieczorek thanks for reporting this problem!

I tested your code, and there seems to be an OS-dependency:

Windows: For me, the problem does not occur, neither using GMT 6.4 nor with GMT 6.5
Linux: I can reproduce this issue already with GMT 6.4
Mac OS: Maybe somebody using Mac OS can check this (:

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()

Linux	Windows

MarkWieczorek commented 3 months ago

Same problem on my intel macOS machine.

seisman commented 3 months ago

I can reproduce your issue. If I change projection='W0/10c' to projection='W10c', it works.

I think it's an upstream bug. Will see if I can fix it later.

yvonnefroehlich commented 3 months ago

I can reproduce your issue. If I change projection='W0/10c' to projection='W10c', it works.

projection="W10c", region="g" or projection="W180/10c", region="d" work (please see the upper row of the figure in my post https://github.com/GenericMappingTools/pygmt/issues/3331#issuecomment-2231416218). Using another central longitude fails (under Linux and Mac OS).

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()

longitude 0°-360°	longitude -180°-180°

seisman commented 2 months ago

The error is caused by an upstream bug, which is fixed in https://github.com/GenericMappingTools/gmt/pull/8554.

seisman commented 2 months ago

The upstream PR https://github.com/GenericMappingTools/gmt/pull/8554 has been merged. With the GMT dev version installed, the following test shows now it works well for any centering longtiude:

import numpy as np
import xarray as xr
import pygmt

# grid that includes redundant 360 E
grid = np.arange(0, 181 * 361).reshape((181, 361))
lat = np.arange(90, -91, -1)
lon = np.arange(0, 361, 1)
grid *= lon

da = xr.DataArray(grid, coords=[('lat', lat), ('lon', lon)])
da.gmt.registration = 0
da.gmt.gtype = 1

# grid that does not include redundant 360 E
grid2 = np.arange(0, 181 * 360).reshape((181, 360))
lat2 = np.arange(90, -91, -1)
lon2 = np.arange(0, 360, 1)
grid2 *= lon2

da2 = xr.DataArray(grid2, coords=[('lat', lat2), ('lon', lon2)])
da2.gmt.registration = 0
da2.gmt.gtype = 1

fig = pygmt.Figure()
for lon in np.arange(0.5, 360, 30):
    fig.grdimage(da, projection=f'W{lon}/10c', region='g', frame=f"+t{lon:.2f} with redundant 360")
    fig.shift_origin(xshift="w+2c")
    fig.grdimage(da2, projection=f'W{lon}/10c', region='g', frame=f"+t{lon:.2f} without redundant 360")
    fig.shift_origin(xshift="-w-2c", yshift="h+2c")
fig.show()

A similar test is added in #3358.

Please note that, if you run the example in https://github.com/GenericMappingTools/pygmt/issues/3331#issue-2411379687, 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.

seisman commented 2 weeks ago

Here is a more comprehensive test for 4 different grids:

da1: Global grid in the longitude [-180, 180] range with redundant longitude at 180/-180
da2: Global grid in the longitude [0, 360] range with redundant longitude at 0/360
da3: Global grid in the longitude [-180, 180] range without redundant longitude at 180/-180
da4: Global grid in the longitude [0, 360] range without redundant longitude at 0/360

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()

The output is: map

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.

GenericMappingTools / pygmt