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

MarkWieczorek · 2024-07-16T14:56:09Z

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 · 2024-07-16T17:06:01Z

@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 · 2024-07-16T19:35:06Z

Same problem on my intel macOS machine.

seisman · 2024-07-17T05:19:59Z

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 · 2024-07-17T17:29:52Z

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 #3331 (comment)). 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 · 2024-07-25T15:05:09Z

The error is caused by an upstream bug, which is fixed in GenericMappingTools/gmt#8554.

seisman · 2024-07-25T15:53:34Z

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:

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

for lon in np.random.rand(30) * 360:
    fig = pygmt.Figure()
    fig.grdimage(da, projection=f'W{lon}/10c', region='g', frame=f"+t{lon:.2f}")
    fig.shift_origin(xshift="w+2c")
    fig.grdimage(da2, projection=f'W{lon}/10c', region='g', frame=f"+t{lon:.2f}")
    fig.show()

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.

MarkWieczorek added the bug Something isn't working label Jul 16, 2024

MarkWieczorek mentioned this issue Jul 16, 2024

plotgmt does not work for grids with extend=False SHTOOLS/SHTOOLS#489

Closed

seisman linked a pull request Jul 25, 2024 that will close this issue

Add a test for checking xr.DataArray grid with Hammer projection #3358

Draft

seisman modified the milestones: 0.13.0, 0.14.0 Jul 25, 2024

seisman removed this from the 0.14.0 milestone Sep 5, 2024

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

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

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

MarkWieczorek commented Jul 16, 2024 •

edited

Loading

yvonnefroehlich commented Jul 16, 2024 •

edited

Loading

MarkWieczorek commented Jul 16, 2024

seisman commented Jul 17, 2024

yvonnefroehlich commented Jul 17, 2024 •

edited

Loading

seisman commented Jul 25, 2024

seisman commented Jul 25, 2024

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

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

Comments

MarkWieczorek commented Jul 16, 2024 • edited Loading

Description of the problem

Minimal Complete Verifiable Example

Full error message

System information

yvonnefroehlich commented Jul 16, 2024 • edited Loading

MarkWieczorek commented Jul 16, 2024

seisman commented Jul 17, 2024

yvonnefroehlich commented Jul 17, 2024 • edited Loading

seisman commented Jul 25, 2024

seisman commented Jul 25, 2024

MarkWieczorek commented Jul 16, 2024 •

edited

Loading

yvonnefroehlich commented Jul 16, 2024 •

edited

Loading

yvonnefroehlich commented Jul 17, 2024 •

edited

Loading