fixed missing crs var in saved network (#719)

* fixed missing crs var in saved network

* simplified calcdist_fromlatlon in polyline_mapslice

* added unittest for polyline_mapslice

* updated whatsnew

dfm_tools/get_nc.py

@@ -153,7 +153,7 @@ def get_xzcoords_onintersection(uds, face_index, crs_dist_starts, crs_dist_stops
     return xr_crs_ugrid
 
 
-def polyline_mapslice(uds:xu.UgridDataset, line_array:np.array, calcdist_fromlatlon:bool = None) -> xu.UgridDataset:
+def polyline_mapslice(uds:xu.UgridDataset, line_array:np.array) -> xu.UgridDataset:
     """
     Slice trough mapdata, combine: intersect_edges_withsort, calculation of distances and conversion to ugrid dataset.
 
@@ -184,15 +184,7 @@ def polyline_mapslice(uds:xu.UgridDataset, line_array:np.array, calcdist_fromlat
     if len(edge_index) == 0:
         raise ValueError('polyline does not cross mapdata')
 
-    #auto determine if cartesian/sperical distance should be computed
-    if calcdist_fromlatlon is None:
-        if hasattr(uds,'projected_coordinate_system'):
-            calcdist_fromlatlon = False
-        elif hasattr(uds,'wgs84'):
-            calcdist_fromlatlon = True
-        else:
-            raise KeyError('To auto determine calcdist_fromlatlon, a variable "projected_coordinate_system" or "wgs84" is required, please provide calcdist_fromlatlon=True/False yourself.')
-    if calcdist_fromlatlon:
+    if uds.grids[0].is_geographic:
         calc_dist = calc_dist_haversine
     else:
         calc_dist = calc_dist_pythagoras

dfm_tools/meshkernel_helpers.py

@@ -138,6 +138,14 @@ def geographic_to_meshkernel_projection(is_geographic:bool) -> meshkernel.Projec
     return projection
 
 
+def meshkernel_is_geographic(mk):
+    if mk.get_projection()==meshkernel.ProjectionType.CARTESIAN:
+        is_geographic = False
+    else:
+        is_geographic = True
+    return is_geographic
+
+
 def meshkernel_to_UgridDataset(mk:meshkernel.MeshKernel, crs:(int,str) = None) -> xu.UgridDataset:
     """
     Convert a meshkernel object to a UgridDataset, including a variable with the crs (used by dflowfm to distinguish spherical/cartesian networks).
@@ -161,6 +169,12 @@ def meshkernel_to_UgridDataset(mk:meshkernel.MeshKernel, crs:(int,str) = None) -
 
     mesh2d_grid = mk.mesh2d_get()
 
+    #check if both crs and grid are geograpic or not
+    #TODO: do this in xugrid: https://github.com/Deltares/xugrid/issues/188
+    grid_is_geographic = meshkernel_is_geographic(mk)
+    if crs_is_geographic != grid_is_geographic:
+        raise ValueError(f"crs has is_geographic={crs_is_geographic} and grid has is_geographic={grid_is_geographic}. This is not allowed.")
+
     # TODO: below is not correctly handled by xugrid yet, projected=False does not give is_geographic=True
     # related issue is https://github.com/Deltares/dfm_tools/issues/686
     xu_grid = xu.Ugrid2d.from_meshkernel(mesh2d_grid, projected= not crs_is_geographic, crs=crs)
@@ -187,7 +201,7 @@ def meshkernel_to_UgridDataset(mk:meshkernel.MeshKernel, crs:(int,str) = None) -
     # add crs including attrs
     if crs is not None:
         xu_grid_uds.ugrid.set_crs(crs)
-
+        uds_add_crs_attrs(xu_grid_uds)
     return xu_grid_uds
 
 
@@ -222,6 +236,7 @@ def uds_add_crs_attrs(uds:(xu.UgridDataset,xr.Dataset)):
     None.
 
     """
+    # TODO: upon crs conversion and to_netcdf(), the netcdf file could get two crs vars, catch this
 
     grid_is_geographic = uds.grid.is_geographic
 
@@ -246,7 +261,9 @@ def uds_add_crs_attrs(uds:(xu.UgridDataset,xr.Dataset)):
     if grid_is_geographic != crs_is_geographic:
         raise ValueError(f"`grid_is_geographic` mismatch between provided grid (is_geographic={grid_is_geographic}) and provided crs ({crs}, is_geographic={crs_is_geographic})")
 
-    #TODO: consider always using the same crs_varn, align with xugrid
+    # TODO: consider always using the same crs_varn, align with xugrid
+    # QGIS also does not recognize epsg anymore when renaming variable to `crs` 
+    # or something else (`wgs84` and `projected_coordinate_system` both do work)
     if grid_is_geographic:
         crs_varn = 'wgs84'
     else:

dfm_tools/xugrid_helpers.py

@@ -160,6 +160,10 @@ def remove_nan_fillvalue_attrs(ds : (xr.Dataset, xu.UgridDataset)):
 
 
 def uds_auto_set_crs(uds : xu.UgridDataset):
+    # FM-mapfiles contain wgs84/projected_coordinate_system variables with epsg attr, xugrid has .crs property
+    # TODO: parse+set crs in xugrid instead: https://github.com/Deltares/xugrid/issues/42
+    # also adjusting projected_coordinate_system/wgs84 when using set_crs/to_crs
+
     uds_epsg = uds.filter_by_attrs(epsg=lambda v: v is not None)
     if len(uds_epsg.data_vars) != 1:
         return
@@ -208,7 +212,6 @@ def open_partitioned_dataset(file_nc, decode_fillvals=False, remove_edges=True,
         - MWRA 3D 20 partitions 2551 timesteps:  74.4/ 3.4 sec (decode_times=False:  79.0 sec)
     
     """
-    #TODO: FM-mapfiles contain wgs84/projected_coordinate_system variables. xugrid has .crs property, projected_coordinate_system/wgs84 should be updated to be crs so it will be automatically handled? >> make dflowfm issue (and https://github.com/Deltares/xugrid/issues/42)
     #TODO: add support for multiple grids via keyword? https://github.com/Deltares/dfm_tools/issues/497
     #TODO: speed up open_dataset https://github.com/Deltares/dfm_tools/issues/225 (also remove_ghost)
 

docs/whats-new.md

@@ -7,6 +7,7 @@
 
 ### Fix
 - prevent ValueError upon concatenation of only emtpy coastlines geodataframes by [@veenstrajelmer](https://github.com/veenstrajelmer) in [#717](https://github.com/Deltares/dfm_tools/pull/717)
+- added dropped crs variable in network file by [@veenstrajelmer](https://github.com/veenstrajelmer) in [#719](https://github.com/Deltares/dfm_tools/pull/719)
 
 
 ## 0.18.0 (2023-12-08)

tests/test_get_nc.py

@@ -0,0 +1,27 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Dec 13 12:17:41 2023
+
+@author: veenstra
+"""
+
+
+import pytest
+import numpy as np
+import dfm_tools as dfmt
+
+
+@pytest.mark.unittest
+def test_polyline_mapslice():
+    uds = dfmt.data.fm_curvedbend_map()
+    timestep = 72
+    line_array = np.array([[ 104.15421399, 2042.7077107 ],
+                            [2913.47878063, 2102.48057382]])
+
+    uds_crs = dfmt.polyline_mapslice(uds.isel(time=timestep), line_array)
+    assert len(uds_crs.grid.node_x) == 720
+    assert uds_crs.grid.face_node_connectivity.shape == (180, 4)
+    assert np.isclose(uds_crs.grid.node_x.min(), 484.02691589067194)
+    assert np.isclose(uds_crs.grid.node_x.max(), 1737.0864257281437)
+    assert np.isclose(uds_crs.grid.node_y.min(), -5)
+    assert np.isclose(uds_crs.grid.node_y.max(), 0.9261683648147339)

tests/test_meshkernel_helpers.py

@@ -154,6 +154,19 @@ def test_geographic_to_meshkernel_projection():
     assert cartesian == cartesian_mk
 
 
+@pytest.mark.systemtest
+def test_meshkernel_to_UgridDataset_geographic_mismatch():
+    """
+    """
+    projection = meshkernel.ProjectionType.CARTESIAN
+    crs = 'EPSG:4326'
+    mk = meshkernel.MeshKernel(projection=projection)
+    try:
+        _ = dfmt.meshkernel_to_UgridDataset(mk=mk, crs=crs)
+    except ValueError as e:
+        assert "crs has is_geographic=True and grid has is_geographic=False" in str(e)
+
+
 @pytest.mark.systemtest
 def test_meshkernel_to_UgridDataset():
     """
@@ -204,6 +217,7 @@ def test_meshkernel_to_UgridDataset():
     assert 0 not in ds_out.mesh2d_face_nodes.to_numpy()
     assert ds_out.mesh2d_face_nodes.to_numpy().min() == -1
     assert ds_out.mesh2d_face_nodes.to_numpy().max() == 626
+    assert "wgs84" in ds_out.variables
 
 
 @pytest.mark.unittest