PyGMT not recognizing numpy ndarray

Q&A PyGMT Q&A
1

I am reading in a raster file using rioxarray. However, when I try to plot it using grdimage, I get the error:
pygmt.exceptions.GMTInvalidInput: Unrecognized data type for grid: <class 'numpy.ndarray'>. How can I fix it? Here is my code:

grid = rioxarray.open_rasterio(‘v9.tif’).drop(‘band’)
fig = pygmt.Figure()
fig.grdimage(grid=grid.data, frame=True)
fig.coast(water=“skyblue”)

2

numpy.ndarray is an accepted input type for functions/methods that expect tabular data. For raster data, the accepted input type is either an xarray.DataArray or a file path (see the grid parameter at https://www.pygmt.org/latest/api/generated/pygmt.Figure.grdimage.html#pygmt.Figure.grdimage).

If you are trying to plot an image rather than a 2D grid, you currently need to provide a file path. The PyGMT team is working on improved support for images loaded as objects in Python (e.g., https://github.com/GenericMappingTools/pygmt/issues/1555).

3

Thank you @maxrjones, this is very helpful! I am trying to plot a raster, but I would prefer to to use images loaded as objects.

4

You can actually plot a GeoTIff directly using fig.grdimage(grid="v9.tif", frame=True), see if it works for your dataset.

5

Thank you, that direct usage of tif file does work great! However, my work case typically involves a raster that has been already been read in and maybe some values altered. I am trying a workaround where I write the raster object into a temporary file on disk using the tempfile library in Python. However, I am have not been able to get that working yet. Do you have a workaround solution for raster objects? Thanks again for some excellent work on this package.

6

Well, it depends on what you’re trying to plot I guess. Is there a reason you’re using fig.grdimage(grid=grid.data, ...) instead of just fig.grdimage(grid=grid)?

The thing is, PyGMT/GMT’s grdimage function requires coordinates in order to know how to label the x-axis and y-axis. These coordinates don’t exist in a numpy.ndarray since it is purely an array of numbers.

What you could try is to put a numpy.ndarray into an xarray.DataArray (i.e. have the coordinates included) following Data Structures. Something like this:

import xarray as xr

height, width = your_numpy_array.shape

new_grid = xr.DataArray(
    data=your_numpy_array,
    coords=dict(
        y=np.linspace(start=miny, stop=maxy, num=height),
        x=np.linspace(start=minx, stop=maxx, num=width),
    ),
    dims=("y", "x"),
)

Or, if your data processing produces an array that is the same shape of the input, you could just use:

new_grid = xr.DataArray(
    data=your_numpy_array,
    coords=dict(y=grid.y, x=grid.x),
    dims=("y", "x"),
)

and then you can just do:

fig.grdimage(grid=new_grid, frame=True)
2 Likes
7

Thank you! Your solution works great! I used grid.data because when I use just grid I get the following error:
ValueError: different number of dimensions on data and dims: 3 vs 2

8

Hmm, that error message looks strange. Could you post the output of print(grid) so that we can see what the original dimensions and data shape of the GeoTiff look like? I found a similar issue at https://github.com/pydata/xarray/issues/3437 but it didn’t seem very helpful.

9

Sure! Here’s the output for print(grid)

<xarray.DataArray (band: 1, y: 3600, x: 7200)>
[25920000 values with dtype=int32]
Coordinates:

  • band (band) int32 1
  • x (x) float64 -180.0 -179.9 -179.9 -179.8 … 179.9 179.9 180.0
  • y (y) float64 89.97 89.92 89.88 89.82 … -89.88 -89.93 -89.98
    spatial_ref int32 0
    Attributes:
    RepresentationType: ATHEMATIC
    STATISTICS_MAXIMUM: 10000
    STATISTICS_MEAN: 461.50334637346
    STATISTICS_MINIMUM: 0
    STATISTICS_SKIPFACTORX: 1
    STATISTICS_SKIPFACTORY: 1
    STATISTICS_STDDEV: 1886.2679169263
    _FillValue: -2147483648.0
    scale_factor: 1.0
    add_offset: 0.0
10

Ah ok, you could try using plotting grid.isel(band=1) instead of grid.data perhaps. That should drop the band dimension and you wouldn’t need to do the extra steps I mentioned above.