563 update meshkernel api calls after new release (#621)

* mk api changes

* update to mk projectiontype

* updated testbank unittests and example scripts

* removed deprecated/commented code from example

* update pyproject.toml and requirements.txt with new hydrolib-core and meshkernel versions

* add geographic_to_meshkernel_projection and use in two functions

dfm_tools/meshkernel_helpers.py

@@ -29,13 +29,13 @@
     "interpolate_bndpli",
     ]
 
-def meshkernel_delete_withcoastlines(mk:meshkernel.meshkernel.MeshKernel, res:str='f', min_area:float = 0, crs:(int,str) = None):
+def meshkernel_delete_withcoastlines(mk:meshkernel.MeshKernel, res:str='f', min_area:float = 0, crs:(int,str) = None):
     """
     Wrapper around meshkernel_delete_withgdf, which automatically gets the bbox from the meshkernel object and retrieves the coastlines_gdf.
 
     Parameters
     ----------
-    mk : meshkernel.meshkernel.MeshKernel
+    mk : meshkernel.MeshKernel
         DESCRIPTION.
     res : str, optional
         DESCRIPTION. The default is 'f'.
@@ -59,13 +59,13 @@ def meshkernel_delete_withcoastlines(mk:meshkernel.meshkernel.MeshKernel, res:st
     meshkernel_delete_withgdf(mk, coastlines_gdf)
 
 
-def meshkernel_delete_withshp(mk:meshkernel.meshkernel.MeshKernel, coastlines_shp:str, crs:(int,str) = None):
+def meshkernel_delete_withshp(mk:meshkernel.MeshKernel, coastlines_shp:str, crs:(int,str) = None):
     """
     Delete parts of mesh that are inside the shapefile polygon.
 
     Parameters
     ----------
-    mk : meshkernel.meshkernel.MeshKernel
+    mk : meshkernel.MeshKernel
         DESCRIPTION.
     coastlines_shp : str
         Path to the shp file.
@@ -88,13 +88,13 @@ def meshkernel_delete_withshp(mk:meshkernel.meshkernel.MeshKernel, coastlines_sh
     meshkernel_delete_withgdf(mk, coastlines_gdb)
 
 
-def meshkernel_delete_withgdf(mk:meshkernel.meshkernel.MeshKernel, coastlines_gdf:gpd.GeoDataFrame):
+def meshkernel_delete_withgdf(mk:meshkernel.MeshKernel, coastlines_gdf:gpd.GeoDataFrame):
     """
     Delete parts of mesh that are inside the polygons/Linestrings in a GeoDataFrame.
 
     Parameters
     ----------
-    mk : meshkernel.meshkernel.MeshKernel
+    mk : meshkernel.MeshKernel
         DESCRIPTION.
     coastlines_gdf : gpd.GeoDataFrame
         DESCRIPTION.
@@ -112,17 +112,17 @@ def meshkernel_delete_withgdf(mk:meshkernel.meshkernel.MeshKernel, coastlines_gd
 
         delete_pol_geom = meshkernel.GeometryList(x_coordinates=xx, y_coordinates=yy) #TODO: .copy()/to_numpy() makes the array contiguous in memory, which is necessary for meshkernel.mesh2d_delete()
         mk.mesh2d_delete(geometry_list=delete_pol_geom, 
-                         delete_option=meshkernel.DeleteMeshOption(2), #ALL_COMPLETE_FACES/2: Delete all faces of which the complete face is inside the polygon
+                         delete_option=meshkernel.DeleteMeshOption.INSIDE_NOT_INTERSECTED,
                          invert_deletion=False)
 
 
-def meshkernel_check_geographic(mk:meshkernel.meshkernel.MeshKernel) -> bool:
+def meshkernel_check_geographic(mk:meshkernel.MeshKernel) -> bool:
     """
     Get projection from meshkernel instance
 
     Parameters
     ----------
-    mk : meshkernel.meshkernel.MeshKernel
+    mk : meshkernel.MeshKernel
         DESCRIPTION.
 
     Returns
@@ -132,13 +132,35 @@ def meshkernel_check_geographic(mk:meshkernel.meshkernel.MeshKernel) -> bool:
 
     """
 
-    if mk.get_projection()==meshkernel.ProjectionType.SPHERICAL:
-        is_geographic = True
-    else:
+    if mk.get_projection()==meshkernel.ProjectionType.CARTESIAN:
         is_geographic = False
+    else:
+        is_geographic = True
     return is_geographic
 
 
+def geographic_to_meshkernel_projection(is_geographic:bool) -> meshkernel.ProjectionType:
+    """
+    converts is_geographic boolean to meshkernel.ProjectionType (SPHERICAL OR CARTESIAN)
+
+    Parameters
+    ----------
+    is_geographic : bool
+        DESCRIPTION.
+
+    Returns
+    -------
+    projection : TYPE
+        DESCRIPTION.
+
+    """
+    if is_geographic:
+        projection = meshkernel.ProjectionType.SPHERICAL
+    else:
+        projection = meshkernel.ProjectionType.CARTESIAN
+    return projection
+
+
 def meshkernel_to_UgridDataset(mk:meshkernel.MeshKernel, crs:(int,str) = None, remove_noncontiguous:bool = False) -> xu.UgridDataset:
     """
     Convert a meshkernel object to a UgridDataset, including a variable with the crs (used by dflowfm to distinguish spherical/cartesian networks).
@@ -268,6 +290,8 @@ def make_basegrid(lon_min,lon_max,lat_min,lat_max,dx,dy,angle=0,is_geographic=Tr
     """
     empty docstring
     """
+    projection = geographic_to_meshkernel_projection(is_geographic)
+
     # create base grid
     make_grid_parameters = meshkernel.MakeGridParameters(angle=angle,
                                                          origin_x=lon_min,
@@ -277,8 +301,8 @@ def make_basegrid(lon_min,lon_max,lat_min,lat_max,dx,dy,angle=0,is_geographic=Tr
                                                          block_size_x=dx,
                                                          block_size_y=dy)
 
-    mk = meshkernel.MeshKernel(is_geographic=is_geographic)
-    mk.curvilinear_make_uniform_on_extension(make_grid_parameters)
+    mk = meshkernel.MeshKernel(projection=projection)
+    mk.curvilinear_compute_rectangular_grid_on_extension(make_grid_parameters)
     mk.curvilinear_convert_to_mesh2d() #convert to ugrid/mesh2d
 
     return mk
@@ -309,14 +333,14 @@ def refine_basegrid(mk, data_bathy_sel, min_edge_size):
     return mk
 
 
-def generate_bndpli_cutland(mk:meshkernel.meshkernel.MeshKernel, res:str='f', min_area:float = 0, crs:(int,str) = None, buffer:float = 0):
+def generate_bndpli_cutland(mk:meshkernel.MeshKernel, res:str='f', min_area:float = 0, crs:(int,str) = None, buffer:float = 0):
     """
     Generate a boundary polyline from the meshkernel object and cut away the landward part.
     Be sure to do this on the base/refined grid, not on the grid where the landward cells were already cut.
     
     Parameters
     ----------
-    mk : meshkernel.meshkernel.MeshKernel
+    mk : meshkernel.MeshKernel
         DESCRIPTION.
     res : str, optional
         DESCRIPTION. The default is 'f'.

dfm_tools/xugrid_helpers.py

@@ -520,25 +520,22 @@ def add_network_cellinfo(uds:xu.UgridDataset):
     """
     #check if indeed 1D grid object
     assert isinstance(uds.grid, xu.Ugrid1d)
-
-    # simple approach, but results in cartesian grid
-    # mk.mesh2d_set(input_mesh2d)
-    # uds_ugrid2d = xu.Ugrid2d.from_meshkernel(mk.mesh2d_get())
-    #TODO: is_geographic=False is currently hardcoded: https://github.com/Deltares/xugrid/issues/128
-
+
     # derive meshkernel from grid with Mesh1d
     mk1 = uds.grid.meshkernel
     mk_mesh1d = mk1.mesh1d_get()
 
     # use Mesh1d nodes/edgenodes info for generation of meshkernel with Mesh2d
     is_geographic = get_uds_isgeographic(uds)
+    from dfm_tools.meshkernel_helpers import geographic_to_meshkernel_projection
+    projection = geographic_to_meshkernel_projection(is_geographic)
     mk_mesh2d = meshkernel.Mesh2d(mk_mesh1d.node_x, mk_mesh1d.node_y, mk_mesh1d.edge_nodes)
-    mk2 = meshkernel.MeshKernel(is_geographic=is_geographic)
+    mk2 = meshkernel.MeshKernel(projection=projection)
     mk2.mesh2d_set(mk_mesh2d)
     mesh2d_grid = mk2.mesh2d_get() #this is a more populated version of mk_mesh2d, needed for xugrid
     #TODO: we have to supply is_geographic twice, necessary?
     # also "projected" is opposite of "is_geographic" according to the docstring
-    xu_grid = xu.Ugrid2d.from_meshkernel(mesh2d_grid, projected= not is_geographic)
+    xu_grid = xu.Ugrid2d.from_meshkernel(mesh2d_grid, projected = not is_geographic)
 
     # convert uds.obj (non-grid vars from dataset) to new xugrid standards
     rename_dims_dict = {uds.grid.node_dimension:xu_grid.node_dimension,

pyproject.toml

@@ -49,10 +49,10 @@ dependencies = [
 	"pydap>=3.4.0",
 	#pooch<1.1.0 do not have attribute retrieve
 	"pooch>=1.1.0",
-	#hydrolib-core<0.5.1 does not yet contain 3D z-sigma mdu settings and many other essentials
-	"hydrolib-core>=0.5.1",
-	#meshkernel<2.1.0 does not yet contain mesh2d_refine_based_on_gridded_samples and other essentials
-	"meshkernel>=2.1.0",
+	#hydrolib-core<0.6.0 does not work with meshkernel>=3.0.0
+	"hydrolib-core>=0.6.0",
+	#meshkernel<3.0.0 api is fundamentally different
+	"meshkernel>=3.0.0",
 ]
 classifiers = [
 	"Development Status :: 4 - Beta",

requirements.txt

@@ -14,5 +14,5 @@ xugrid>=0.7.0
 cdsapi>=0.6.1
 pydap>=3.4.0
 pooch>=1.1.0
-hydrolib-core>=0.5.1 
-meshkernel>=2.1.0
+hydrolib-core>=0.6.0
+meshkernel>=3.0.0

tests/examples_workinprogress/workinprogress_meshkernel_creategrid.py

@@ -10,17 +10,15 @@
 import xarray as xr
 import matplotlib.pyplot as plt
 plt.close('all')
-import numpy as np
 import contextily as ctx
 import dfm_tools as dfmt
-
-#TODO: maybe use make_basegrid and refine_basegrid functions from dfmt.meshkernel_helpers
+import datetime as dt
 
 #general settings
 lon_min,lon_max = -6,2
 lat_min,lat_max = 48.5,51.2
 lon_res,lat_res = 0.2,0.2
-is_geographic = True #True for spherical grids
+projection = meshkernel.ProjectionType.SPHERICAL #True for spherical grids
 
 figsize = (10,4)
 crs = 'EPSG:4326'
@@ -33,6 +31,7 @@
 
 """
 
+#TODO: maybe use make_basegrid and refine_basegrid functions from dfmt.meshkernel_helpers
 # Create an instance of MakeGridParameters and set the values
 make_grid_parameters = meshkernel.MakeGridParameters(angle=0.0, #TODO: does non-zero result in an orthogonal spherical grid?
                                                      origin_x=lon_min,
@@ -42,19 +41,10 @@
                                                      block_size_x=lon_res,
                                                      block_size_y=lat_res)
 
-grid_in_pol = False
-# If a polygon is provided it will be used in the generation of the curvilinear grid. The polygon must be closed
-if grid_in_pol: #can be used instead of origin_x/origin_y and num_x/num_y
-    pol_x = np.array([-6,-4,0,-6], dtype=np.double)
-    pol_y = np.array([48,51,49.5,48], dtype=np.double)
-    geometry_list = meshkernel.GeometryList(pol_x, pol_y)
-else:
-    geometry_list = None
-
-
-mk = meshkernel.MeshKernel(is_geographic=is_geographic)
-mk.curvilinear_make_uniform_on_extension(make_grid_parameters) #TODO: geometry_list is possible with curvilinear_make_uniform, but not for this function: https://github.com/Deltares/MeshKernelPy/issues/76
+mk = meshkernel.MeshKernel(projection=projection)
+mk.curvilinear_compute_rectangular_grid_on_extension(make_grid_parameters)
 mk.curvilinear_convert_to_mesh2d() #convert to ugrid/mesh2d
+
 mesh2d_basegrid = mk.mesh2d_get() #in case of curvi grid: mk.curvilinear_convert_to_mesh2d()
 fig, ax = plt.subplots(figsize=figsize)
 mesh2d_basegrid.plot_edges(ax,linewidth=0.8)
@@ -79,7 +69,6 @@
 values_np = data_bathy_sel.elevation.to_numpy().flatten().astype('float') #TODO: astype to avoid "TypeError: incompatible types, c_short_Array_74880 instance instead of LP_c_double instance"
 gridded_samples = meshkernel.GriddedSamples(x_coordinates=lon_np,y_coordinates=lat_np,values=values_np) #TODO: does not result in refinement
 
-
 #refinement
 mesh_refinement_parameters = meshkernel.MeshRefinementParameters(#refine_intersected=False, #TODO: what does this do?
                                                                  #use_mass_center_when_refining=False, #TODO: what does this do?
@@ -106,30 +95,17 @@
 delete (landward) part of grid with polygon and plot result
 """
 
-#line, = ax.plot([], [],'o-') # empty line
-# #linebuilder = dfmt.LineBuilder(line) #TODO: this makes it possible to interactively click a line in the bedlevel figure. Use linebuilder.line_array as alternative line_array
-# delete_pol = np.array([[ 1.91741935, 49.76580645],
-#                         [ 0.20387097, 49.9       ],
-#                         [-0.25032258, 48.71290323],
-#                         [ 1.92774194, 48.59935484]])
 print('deleting landward grid with coastlines (takes a while)...')
+dtstart = dt.datetime.now()
 dfmt.meshkernel_delete_withcoastlines(mk=mk, res='h', min_area=100)
-print('done')
-
+print(f'done in {(dt.datetime.now()-dtstart).total_seconds():.2f} sec')
 
 mesh2d_noland = mk.mesh2d_get()
 fig, ax = plt.subplots(figsize=figsize)
 mesh2d_noland.plot_edges(ax,linewidth=0.8)
-# xlim,ylim = ax.get_xlim(),ax.get_ylim() #get x/ylims before ldb plotting changes it
-# for iP, pol_del in enumerate(pol_ldb_list):
-#     ax.plot(pol_del['x'],pol_del['y'],'-r')
-# ax.set_xlim(xlim)
-# ax.set_ylim(ylim)
 dfmt.plot_coastlines(ax=ax, res='h', min_area=100)
 ctx.add_basemap(ax=ax, crs=crs, attribution=False)
 
-
-
 """
 convert meshkernel grid to xugrid, interp bathymetry, plot and save to *_net.nc
 """
@@ -138,6 +114,7 @@
 
 fig, ax = plt.subplots(figsize=figsize)
 xu_grid_uds.grid.plot(ax=ax,linewidth=0.8) #TODO: maybe make uds instead of ds (but then bathy interpolation goes wrong)
+dfmt.plot_coastlines(ax=ax, res='h', min_area=100)
 ctx.add_basemap(ax=ax, crs=crs, attribution=False)
 
 #interp bathy
@@ -153,7 +130,5 @@
 netfile = 'englishchannel_net.nc'
 xu_grid_uds.ugrid.to_netcdf(netfile)
 
-#TODO: update https://github.com/Deltares/dfm_tools/issues/217
-
 #TODO: small flow links in resulting grid