cell center vs. cell corner confusion when using Local Coordinates

[Snowthe]

[This might be related to https://github.com//issues/540 or the fix of it]

Describe the bug
When using dfs2/3 files with local coordinates, one defines the corners coordinates of the origin cell in its properties - map projection coordinates - Easting and Northing. This is relfected correctly when viewing a dfs2/3 file in the MIKE Zero Result viewer, where one can see the correct corner coordinates.
When I load the same file with mikeio.read(), and inspect the mikeio.Dataset.geometry, it shows the same corner coordinates - as expected.
However, as soon as I try to plot the array (or write to tif via rioxarray), I can see that it is shifted by half a cell in both directions - i.e. the corner coordinates are wrongly interpreted as center coordinates.

To Reproduce

fp = 'test.dfs2'
ds = mikeio.read(fp)
ds
ds['item'].plot()

Is this related to some inconsistency between a mikeio Dataset/DataArray and the respective xarray data types?

System information:

• Python version 3.8
• MIKE IO version 1.6.1

[ecomodeller]

It seems like we are MIKE IO is always considering the coordinates to be cell centers irrespective of CRS.

According to:
MIKE FAQ on this topic https://faq.dhigroup.com/default.asp?module=MIKE+Zero&ID=284
it seems like we should treat "NON-UTM" differently.

@jsmariegaard do we have any documentation other than the FAQ where this convention is clearly stated?

[Snowthe]

Thanks for the quick follow-up!
With "we are always considering the coordinates to be cell centres", you mean in mikeio? Which might differ from MIKE Zero as such?

[jsmariegaard]

Yes, I know. For me it is a known issue - but that should of course be communicated somehow and even better fixed 🙄

It was not easy to make an elegant solution to this inconsistency.

[Snowthe]

...so in case I would like to use a workaround in the meantime - what would be a good way?
I am not aware of a possiblity of directly setting an attribute (such as Grid2D.origin) - i.e. I would create a new Grid2D/3D geometry with modified origin, and then "reassemble" the dataset as shown in this section: https://dhi.github.io/getting-started-with-mikeio/dfs2.html#writing-data ?

[jsmariegaard]

We are doing something similar when we convert a Grid2D to a mesh - see

mikeio/mikeio/spatial/_grid_geometry.py

Line 921 in 66beb8d

def to_geometryFM(self, *, z=None, west=2, east=4, north=5, south=3):

I don't know of this or any of the helper functions here could help?

[Snowthe]

This issue seems to be unresolved as of today - and I really would appreciate it a fix.

1. As long as you fully remain in mikeio (e.g. open a dfs2 file via mikeio, do some calculations on the mikeio dataset, save the output to another dfs2), all is fine.
2. That, however, does not apply to the plotting method in mikeio; that one also confuses center and corner coordinates.
3. However, as soon as you move into e.g. something with xarray (ds.to_xarray() to perform some more advanced calculations), or save it to a geotiff or similar, you need to correct for the center/corner confusion.

I have written a helper function to correct for it in case 3; however you may not apply it in case 1, so it easily gets confusing.

[ecomodeller]

@Snowthe Is this the behaviour you need?

[Snowthe]

@Snowthe Is this the behaviour you need?

Yes, awesome, this looks like the solution. (As of the current version, this results in -0.25, -0.25 as corner coordinates instead of 0,0)

I assume that this not only applies to da.plot(), but is properly passed on in case one e.g. uses xr_da = da.to_xarray() (and then work on with the resulting xarray, plot it etc)

[ecomodeller]

@Snowthe Would you be so kind to verify that this is the solution you need?

You can install the branch like this:

pip install https://github.com/DHI/mikeio/archive/non-utm.zip

[Snowthe]

@Snowthe Would you be so kind to verify that this is the solution you need?

You can install the branch like this:

pip install https://github.com/DHI/mikeio/archive/non-utm.zip

Thanks a lot. To me, it seems as if this gives the desired behaviour. I checked

import mikeio
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import rasterio
import rioxarray

da = mikeio.DataArray(data=np.array([[1,2],[3,4]]),
                       geometry=mikeio.Grid2D(dx=0.5, projection='NON-UTM', nx=2, ny=2))

da.geometry.x #now 0.25, 0.75
plt.figure(); da.plot() #now lower left corner as 0,0 (instead -0.25, -0.25)

xr_da = da.to_xarray()
xr_da.x #now 0.25, 0.75
plt.figure(); xr_da.plot() #now lower left corner as 0,0 (instead -0.25, -0.25)

xr_da.rio.to_raster('testtif2.tif') #now lower left corner as 0,0 (instead -0.25, -0.25)

[Snowthe]

I fear I have to come back to this - I just realized that loading dfs2 files seems to work fine with the 2.0.dev0 version you provided me. Loading dfs3 files, however, still seems to show the same error in some cases:

ds = mikeio.read(dfs3_fp)
Seems to work fine. However, when I specify a specific layer via
ds = mikeio.read(dfs3_fp, layer=1)
It seems to go wrong as before

[Snowthe]

Thanks for re-opening.
As an update: ds = mikeio.read(dfs3_fp, layer=1) seems to be correcting wrongly (i.e. in the wrong direction), instead of not correcting at all (=old behaviour) for the NON-UTM center/corner confusion.