Issue with DEM()

[filippogiaroli]

Hi!

I am trying to figure out the different uses of DEM() since I have a netcdf that is not working. So I tried both the string and array options with an example that I know it works.

I have a netcdf, that if I load and plot with xarray:

plt.figure()
plt.contourf(ds.lon.data, ds.lat.data, ds.z.data)
plt.colorbar()

If I then use DEM() it all works correctly:

dem_etopo_built_nc = f"/home/filippog/Mesh/dem_etopo_built_{name}.nc"
dem_etopo_fine = om.DEM(dem_etopo_built_nc)
dem_etopo_fine.plot()

I wanted to get the same result giving an array as input.
1st ISSUE: you need to explicitly give bbox to DEM() since it fails to create it inside the function if not given (bbox=None) and then it crashes in lines 677-679 of geodata.py

bbox = (8, 10, 44, 45) # Liguria
dem_etopo_fine = om.DEM(ds.z.values, bbox=(8, 10, 44, 45))
dem_etopo_fine.plot()

I get a wrong dem:

After playing with the array I figured out that you need to give it as:
dem_etopo_fine = om.DEM(ds.z.values[::-1, :], bbox=(8, 10, 44, 45))
But cannot figure out why I need to invert the lon?

Any thoughts or ideas?

Thanks a lot!!

[SorooshMani-NOAA]

The main issue when you pass in an array of only z values is that you'd lose all the information about location. So the bbox actually defines locations for you, that's why instead of clipping the original DEM with the given bbox, you actually get the stretched/shrunk version of the same original elevations. The issue with the inversion is related to how data is stored vs plotted when you lose the location info. Your first row of data in the matrix starts from the top. In the plot, your first row of pixels starts from the bottom, and it looks like that's what you see here.

I haven't looked at the DEM code in oceanmesh, I've recently started looking into using this package. So please don't take what I said as given!

[ml14je]

Hi there. Thank you for flagging this! I think I was the one who initially implemented the input of an array/callable function for DEM. It's more or less what was said above in relation to needing to prescribe bbox: by giving an array, you need some information about the positioning. The tuple bbox is required as it defines the positional data. For the input of numpy data, there is also a supposition of how that data is given. It could be that a np.rot90 is needed just after line 676 of geodata.py. I can look into this Wednesday evening. If I understand correctly, you do have a fix for the time-being?

[krober10nd]

From memory, DEM class takes a tif or netcdf file and reads it in as you show.

I would avoid passing arrays with no meta data.

Sometimes the NetCdf/tif file is packed using a different convention than is expected by the software. What I recall is that the origin is either at the top left or bottom left of the DEM.

In any case, this aspect of the software should be improved and standardized.

[krober10nd]

Hi @filippogiaroli I have finally found some time to improve this issue. I will merge the branch referenced above shortly but if you want you can use that branch now. Please let me know if you run into any issues with reading or processing your DEM with that branch if you try it. Thanks.