csdms / pymt

A Python toolkit for running and coupling Earth surface models
https://csdms.colorado.edu/wiki/PyMT
MIT License
51 stars 19 forks source link

Incorrect dimension used for rectilinear grid in _BmiCap grid methods #135

Closed mdpiper closed 3 years ago

mdpiper commented 3 years ago

In the _BmiBap.grid_x method, an array is allocated to hold the x values of a rectilinear grid:

https://github.com/csdms/pymt/blob/79855ce2393c382d3f7b0e22b0d5a55fd582fe21/pymt/framework/bmi_bridge.py#L460-L461

The indexing holds for arrays of rank 2 and shape (y, x), but not for arrays of rank 3 and shape (z, y, x). The same indexing problem occurs in _BmiBap.grid_y and _BmiBap.grid_z.

The clue is in the traceback from a rank 3 variable:

bmi-geotiff
<xarray.Rectilinear>
Dimensions:     (rank: 3)
Dimensions without coordinates: rank
Data variables:
    mesh        int64 0
    node_shape  (rank) int32 1 1623 2043
Traceback (most recent call last):
  File "/Users/mpiper/build/pymt_geotiff/examples/pymt_geotiff-pymt_ex.py", line 15, in <module>
    m.initialize(*args)
  File "/Users/mpiper/anaconda3/envs/pymt/lib/python3.9/site-packages/pymt/framework/bmi_bridge.py", line 324, in initialize
    self._grid[grid_id] = dataset_from_bmi_grid(self, grid_id)
  File "/Users/mpiper/anaconda3/envs/pymt/lib/python3.9/site-packages/pymt/framework/bmi_ugrid.py", line 318, in dataset_from_bmi_grid
    grid = Rectilinear(bmi, grid_id)
  File "/Users/mpiper/anaconda3/envs/pymt/lib/python3.9/site-packages/pymt/framework/bmi_ugrid.py", line 244, in __init__
    nodes = self.get_nodes()
  File "/Users/mpiper/anaconda3/envs/pymt/lib/python3.9/site-packages/pymt/framework/bmi_ugrid.py", line 48, in get_nodes
    data = getattr(self.bmi, "grid_" + dim_name)(self.grid_id)
  File "/Users/mpiper/anaconda3/envs/pymt/lib/python3.9/site-packages/pymt/framework/bmi_bridge.py", line 465, in grid_x
    self.bmi.get_grid_x(grid, out)
  File "/Users/mpiper/anaconda3/envs/pymt/lib/python3.9/site-packages/bmi_geotiff/bmi.py", line 319, in get_grid_x
    x[:] = self._da.x.values
ValueError: could not broadcast input array from shape (2043,) into shape (1623,)

Note the x and y dimensions are swapped.

mdpiper commented 3 years ago

I'm totally to blame for this issue, but I have a nice solution that I'll propose in a PR.