-
Notifications
You must be signed in to change notification settings - Fork 11
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Add function to modelbuilder to generate observations along coastline #578
base: main
Are you sure you want to change the base?
Changes from all commits
4a0af1d
338949f
5e06f83
a394f39
eafd513
File filter
Filter by extension
Conversations
Jump to
Diff view
Diff view
There are no files selected for viewing
Original file line number | Diff line number | Diff line change |
---|---|---|
|
@@ -9,15 +9,139 @@ | |
""" | ||
|
||
import os | ||
from pathlib import Path | ||
import xarray as xr | ||
import pandas as pd | ||
import dfm_tools as dfmt | ||
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 | ||
from hydrolib.core.dimr.models import DIMR, FMComponent, Start | ||
import warnings | ||
|
||
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): | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Concat all in multilinesrring or split all in linestring, so we prevent this if-else loop. I think you can use the private function |
||
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) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. If you provide the netfile, I suggest to also use |
||
net_face_x = [] | ||
net_face_y = [] | ||
net_face_z = [] | ||
for face in ds_net.mesh2d_nFaces: | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Gebruik hier |
||
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}) | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Eerst zvar ophalen, filter op diepte (sel) en alleen die xy coords bewaren |
||
|
||
#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): | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Might be related to |
||
""" | ||
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): | ||
""" | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Same as for xlonylat2xyzcartesian |
||
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 | ||
There was a problem hiding this comment. Choose a reason for hiding this commentThe reason will be displayed to describe this comment to others. Learn more. Buiten functie halen |
||
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 | ||
|
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,14 @@ | ||
from dfm_tools import modelbuilder as mb | ||
|
||
netfile = r'p:\11208614-de-370a\01_models\BES\DFLOWFM\Models\bonaire\bonaire_net.nc' | ||
|
||
interval = 0.05 | ||
resolution = 'h' | ||
threshold_mindepth = 20 | ||
|
||
obs = mb.generate_coastline_obspoints(interval = interval, | ||
resolution = resolution, | ||
threshold_mindepth = threshold_mindepth, | ||
file_nc = netfile) | ||
|
||
|
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Possible to do this on coastlines_gdb already? Would prevent many lines of code