Add function to modelbuilder to generate observations along coastline

Fixes #575

This function is now added, but needs TLC and cleaning up. To do for Tammo

dfm_tools/modelbuilder.py

@@ -15,7 +15,130 @@
 import hydrolib.core.dflowfm as hcdfm
 import datetime as dt
 import glob
+import geopandas as gpd
+import numpy as np
+import pyproj
+import shapely
+from shapely.geometry import LineString, MultiPoint, Point
+from scipy.spatial import KDTree
 
+def generate_coastline_obspoints(lon_min : float, # degrees
+                                 lon_max : float, # degrees
+                                 lat_min : float, # degrees
+                                 lat_max : float, # degrees
+                                 interval : float,# degrees
+                                 resolution : str ,# degrees
+                                 threshold_mindepth : float, # meters
+                                 file_nc : Path, # path to netcdf file
+                                 output_file : Path): # path to output file,):
+
+    #%% create obs file and snap to grid 
+    # regular  grid in coast
+    lat =  np.linspace(lat_min,lat_max,int((lat_max-lat_min)/0.5))
+    lon =  np.linspace(lon_min,lon_max,int((lon_max-lon_min)/0.5))
+
+    # Points along coastline with dfm_tools
+    coastline = dfmt.get_coastlines_gdb(res = resolution,
+                             bbox = [lon_min, lat_min, lon_max, lat_max],
+                             crs = "EPSG:4326")
+    linestrings = [polygon.boundary for polygon in coastline['geometry']]
+    linstrings_cropped = [shapely.ops.clip_by_rect(linestring, lon_min, lat_min, lon_max, lat_max) for linestring in linestrings]
+    shp_clipped = gpd.GeoDataFrame(geometry = gpd.GeoSeries(linstrings_cropped))
+    shp_clipped = shp_clipped[~shp_clipped.is_empty]
+
+    #Build GeoSeries with points along coastline 
+    cpoints = gpd.GeoSeries()
+    for index, line in shp_clipped.iterrows():
+        if 'MULTILINESTRING' in str(line.values):       
+            shapes = str(line.values[0]).strip().split("\n")       
+            gdf = gpd.GeoDataFrame({'geometry': shapes})
+            gdf['geometry'] = gpd.GeoSeries.from_wkt(gdf['geometry'])
+            gdf = gdf.set_geometry('geometry').explode(index_parts=True)
+            for iindex, iline in  gdf.iterrows():
+                distances = np.arange(0, iline.geometry.length, interval)
+                points = MultiPoint([LineString(iline.geometry).interpolate(distance) for distance in distances])
+                gs =gpd.GeoSeries(Point(pnt.x,pnt.y) for pnt in points.geoms)
+                cpoints=pd.concat([cpoints, gs])   #.append(gs)
+        else:
+            distances = np.arange(0, line.geometry.length, interval)
+            points = MultiPoint([LineString(line.geometry).interpolate(distance) for distance in distances])
+            gs =gpd.GeoSeries(Point(pnt.x,pnt.y) for pnt in points.geoms)
+            cpoints=pd.concat([cpoints, gs])# cpoints=cpoints.append(gs)
+    locs_coast = gpd.GeoDataFrame(cpoints, geometry=cpoints.geometry, crs="EPSG:4326")
+
+    # Snap to model points and -5 m depth ## TO DO
+    obs=pd.DataFrame()
+    obs['x'] = locs_coast.geometry.x.values
+    obs['y'] = locs_coast.geometry.y.values
+
+    #%% calculate face center x, y and z
+    print('interpolating cell centers from nodes')
+    ds_net = xr.open_dataset(file_nc)
+    net_face_x = []
+    net_face_y = []
+    net_face_z = []
+    for face in ds_net.mesh2d_nFaces:
+        idx_node = ds_net.mesh2d_face_nodes\
+            .sel(mesh2d_nFaces = face)\
+            .dropna(dim='mesh2d_nMax_face_nodes').data - 1
+        idx_node =  [int(idx_node[i]) for i in range(len(idx_node))]
+        net_face_x_sel = ds_net.mesh2d_node_x.sel(mesh2d_nNodes = idx_node).data.mean()
+        net_face_y_sel = ds_net.mesh2d_node_y.sel(mesh2d_nNodes = idx_node).data.mean()
+        net_face_z_sel = ds_net.mesh2d_node_z.sel(mesh2d_nNodes = idx_node).data.mean()
+        net_face_x.append(net_face_x_sel)
+        net_face_y.append(net_face_y_sel)
+        net_face_z.append(net_face_z_sel)
+
+    #Make dictionary with cell centers
+    net_faces = pd.DataFrame({'x':net_face_x, 'y':net_face_y, 'z':net_face_z})
+
+    #Select only cells with z < threshold_mindepth
+    bool_valid_cells = (net_faces['z']<-threshold_mindepth)
+
+    #creating kdtree with valid cell centers (cartesian coordinates)
+    def xlonylat2xyzcartesian(data):
+        """
+        necessary to calculate cartesian distances, otherwise nearest neigbour can fail.
+        https://stackoverflow.com/questions/45127141/find-the-nearest-point-in-distance-for-all-the-points-in-the-dataset-python
+        """
+        R = 6367
+        phi = np.deg2rad(data['y'])
+        theta = np.deg2rad(data['x'])
+        data = pd.DataFrame()
+        data['x_cart'] = R * np.cos(phi) * np.cos(theta)
+        data['y_cart'] = R * np.cos(phi) * np.sin(theta)
+        data['z_cart'] = R * np.sin(phi)
+        return data
+
+    data_celcenxy_valid = net_faces.loc[bool_valid_cells, ['x', 'y']].reset_index() #pd.DataFrame({'x':net_node_x[bool_valid_cells],'y':net_node_y[bool_valid_cells]})#,'area':data_cellarea[bool_valid_cells]})
+    data_celcenxy_valid_cart = xlonylat2xyzcartesian(data_celcenxy_valid)
+    tree = KDTree(data_celcenxy_valid_cart[['x_cart','y_cart','z_cart']])
+
+    def dist_to_arclength(chord_length):
+        """
+        https://stackoverflow.com/questions/45127141/find-the-nearest-point-in-distance-for-all-the-points-in-the-dataset-python
+        """
+        R = 6367 # earth radius
+        central_angle = 2*np.arcsin(chord_length/(2.0*R)) 
+        arclength = R*central_angle*1000
+        return arclength
+
+    #finding nearest cellcenter-neighbors of each obspoint in file
+    data_obsorg_cart = xlonylat2xyzcartesian(obs)
+    distance_cart, index = tree.query(data_obsorg_cart, k=1)
+    data_celcenxy_validsel = data_celcenxy_valid.loc[index,:]
+
+    # write 
+    obs_snapped=pd.DataFrame()
+    obs_snapped['x'] = data_celcenxy_validsel['x'].values
+    obs_snapped['y'] = data_celcenxy_validsel['y'].values
+    obs_snapped = obs_snapped.drop_duplicates()
+
+    #Create obs file for dflowfm
+    obs_sfincs = hcdfm.XYNModel()
+    points = [hcdfm.XYNPoint(x=x, y=y, n=f'sfincs_bnd_{i:03d}') for i, (x, y) in enumerate(zip(obs_snapped['x'].values, obs_snapped['y'].values))]
+    obs_sfincs.points = points
+    return obs_sfincs
 
 def cmems_nc_to_bc(ext_bnd, list_quantities, tstart, tstop, file_pli, dir_pattern, dir_output, refdate_str):
     #input examples in https://github.com/Deltares/dfm_tools/blob/main/tests/examples/preprocess_interpolate_nc_to_bc.py