Deltares / xugrid

Xarray and unstructured grids
https://deltares.github.io/xugrid/
MIT License
58 stars 8 forks source link

Back and forth conversion xugrid - meshkernel does not work for spherical network #177

Closed Tammo-Zijlker-Deltares closed 8 months ago

Tammo-Zijlker-Deltares commented 9 months ago

Loading a ugriddataset, converting to a meshkernel object (to apply transformations such as deletion with meshkernel) and converting back to a ugrid dataset does not work for spherical networks.

This probably has to do that in conversion from ugrid to meshkernel, the sperical projection information is lost.

Sample code:

# Imports
import dfm_tools as dfmt
import meshkernel as mk

# Load a model grid as UgridDataset

network = dfmt.open_partitioned_dataset(r'p:\11208614-de-370a\01_models\Humber\DFLOWFM\Models\humber_basque_modelbuilder\humber_basque_net.nc') 

#Convert the network to a meshkernel object
mk_object = network.grid.meshkernel

# Delete non-contiguous parts of the grid
mk_object.mesh2d_remove_disconnected_regions()
crs = 'EPSG:4326'

#convert to xugrid, interpolate z-values and write to netcdf
xu_grid_uds = dfmt.meshkernel_to_UgridDataset(mk=mk_object, crs=crs)

The error is:

ValueError: is_geographic mismatch between provided grid (is_geographic=False) and provided crs (EPSG:4326, is_geographic=True)

veenstrajelmer commented 9 months ago

This has to do with the fact that all Ugrid2d.meshkernel instances are initialized as cartesian. We should support spherical networks in the xugrid to meshkernel conversion.

It is also related to https://github.com/Deltares/dfm_tools/issues/661 and https://github.com/Deltares/xugrid/issues/42. If we set the uds crs automatically, we can derive is_geographic within uds.grid.meshkernel and apply cartesian/spherical projection correctly.

Huite commented 9 months ago

Related to: https://github.com/Deltares/xugrid/issues/52

Any code that touches lengths or areas requires an alternative implementation.

It shouldn't be hard to support initializing things, but the next question is whether should put up some guard rails. I.e. not allowing certain operations if the cartesian flag on the geometry is False. Maybe just raise some NotImplementedErrors.

veenstrajelmer commented 8 months ago

With https://github.com/Deltares/xugrid/pull/186 the meshkernel projection is now automatically derived from the uds. It is important to set the crs on the uds however. In the future we will automatically set it in https://github.com/Deltares/xugrid/issues/42 or https://github.com/Deltares/dfm_tools/issues/684

This works:

import dfm_tools as dfmt
crs = 'EPSG:4326'
network = dfmt.open_partitioned_dataset(r'p:\11208614-de-370a\01_models\Humber\DFLOWFM\Models\humber_basque_modelbuilder\humber_basque_net.nc') 
network.ugrid.set_crs(crs)
mk_object = network.grid.meshkernel
xu_grid_uds = dfmt.meshkernel_to_UgridDataset(mk=mk_object, crs=crs)