+
+
+
+
+
+
+
+In [ ]:
+
+
+
+
+
+import glob
+import os
+import io
+import pymap3d
+
+import cartopy.feature as cf
+import geopandas as gpd
+import holoviews as hv
+import hvplot.pandas
+import hvplot.xarray
+import numpy as np
+import pandas as pd
+import shapely
+import thalassa
+import xarray as xr
+
+import searvey
+
+import pyposeidon
+import pyposeidon.meteo as pmeteo
+import pyposeidon.dem as pdem
+import pyposeidon.boundary as pbound
+import pyposeidon.mesh as pmesh
+import pyposeidon.model as pmodel
+from pyposeidon.utils import cast
+
+import pyposeidon.utils.hplot as hplot
+import pyposeidon.utils.pplot as pplot
+
+hv.extension("bokeh")
+
+!mkdir -p data/iceland
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+Geometry¶
+
+
+
+
+
+
+
+
+In [ ]:
+
+
+
+
+
+lon_min = -28.0
+lon_max = -11.1
+lat_min = 62.0
+lat_max = 68.0
+
+bbox = shapely.box(lon_min, lat_min, lon_max, lat_max)
+geometry = dict(lon_min=lon_min, lon_max=lon_max, lat_min=lat_min, lat_max=lat_max)
+
+
+
+
+
+
+
+
+
+Coastlines¶
+
+
+
+
+
+
+
+In [ ]:
+
+
+
+
+
+OSM_FOLDER = "/home/tomsail/work/python/seareport_org/coastlines/raw/osm/land-polygons-complete-4326"
+
+
+
+
+
+
+
+
+In [ ]:
+
+
+
+
+
+coastlines = gpd.read_file(OSM_FOLDER + "/land_polygons.shp", bbox=bbox)
+coastlines
+
+
+
+
+
+
+
+
+Out[Â ]:
+
+
+
+
+
+
+
+
+| + | FID | +geometry | +
|---|---|---|
| 0 | +29026 | +POLYGON ((-22.02066 64.16396, -22.02093 64.163... | +
| 1 | +29034 | +POLYGON ((-21.76014 64.16271, -21.76026 64.162... | +
| 2 | +29035 | +POLYGON ((-21.77483 64.18337, -21.77487 64.183... | +
| 3 | +29062 | +POLYGON ((-21.83118 64.18681, -21.83212 64.187... | +
| 4 | +29065 | +POLYGON ((-21.9012 64.25611, -21.90095 64.2561... | +
| ... | +... | +... | +
| 3637 | +766100 | +POLYGON ((-22.40777 65.09126, -22.40769 65.091... | +
| 3638 | +766101 | +POLYGON ((-22.41448 65.09122, -22.41439 65.091... | +
| 3639 | +769143 | +POLYGON ((-15.21416 64.27749, -15.21428 64.277... | +
| 3640 | +770114 | +POLYGON ((-23.64125 64.73958, -23.64099 64.739... | +
| 3641 | +772163 | +POLYGON ((-19.00579 63.41416, -19.00584 63.414... | +
3642 rows × 2 columns
+
+
+
+
+
+
+
+
+
+Boundaries¶
+
+
+
+
+
+
+
+In [ ]:
+
+
+
+
+
+boundary = pbound.Boundary(geometry=geometry, coastlines=coastlines)
+boundary.contours.head()
+len(boundary.contours)
+
+
+
+
+
+
+
+
+
+
+
+2024-07-15 07:22:41,215 WARNING pyposeidon.boundary __init__:67 coastlines is not a file, trying with geopandas Dataset ++
+
+Out[Â ]:
+
+
+
+
+
+
+
+
+| + | geometry | +tag | +lindex | +nps | +
|---|---|---|---|---|
| 0 | +LINESTRING (-28 68, -11.1 68, -11.1 62, -28 62... | +open | +1 | +4 | +
| 1 | +LINESTRING (-13.85278 65.57077, -13.85307 65.5... | +island | +-1 | +5 | +
| 2 | +LINESTRING (-13.84513 65.56995, -13.845 65.570... | +island | +-2 | +6 | +
| 3 | +LINESTRING (-13.68221 65.5536, -13.68187 65.55... | +island | +-3 | +6 | +
| 4 | +LINESTRING (-13.68278 65.55315, -13.68273 65.5... | +island | +-4 | +5 | +
+
+Out[Â ]:
+
+
+3643+
+
+
+
+
+
+
+In [ ]:
+
+
+
+
+
+boundary.show()
+
+
+
+
+
+
+
+
+Out[Â ]:
+
+
+<Axes: >+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+In [ ]:
+
+
+
+
+
+# boundary.contours.hvplot(geo=True, tiles=True, frame_height=500)
+
+
+
+
+
+
+
+
+
+DEM¶
+
+
+
+
+
+
+
+In [ ]:
+
+
+
+
+
+# url = "https://coastwatch.pfeg.noaa.gov/erddap/griddap/srtm30plus"
+url = "data/iceland.nc"
+dem = pdem.Dem(dem_source=url, **geometry)
+dem.Dataset.load()
+
+
+
+
+
+
+
+
+
+
+
+2024-07-15 07:22:43,862 INFO pyposeidon.tools open_dataset:350 extracting dem from data/iceland.nc + +2024-07-15 07:22:43,863 DEBUG pyposeidon.tools open_dataset:361 Engine: netcdf4 ++
+
+
+
+
+/home/tomsail/miniconda3/envs/pos_test/lib/python3.11/site-packages/pyposeidon/dem.py:131: UserWarning: There are multiple data_vars. Assuming 'elevation' is the first one: ['Band1', 'crs']
+ warnings.warn(f"There are multiple data_vars. Assuming 'elevation' is the first one: {data_vars}")
+
+
+
+
+
+
+2024-07-15 07:22:44,594 INFO pyposeidon.dem dem_:210 dem done + +2024-07-15 07:22:44,595 WARNING pyposeidon.dem __init__:81 coastlines not present, aborting adjusting dem + ++
+
+Out[Â ]:
+
+
+
+
+<xarray.Dataset> Size: 23MB +Dimensions: (latitude: 1444, longitude: 4059) +Coordinates: + * latitude (latitude) float64 12kB 61.99 62.0 62.0 62.01 ... 68.0 68.0 68.01 + * longitude (longitude) float64 32kB -28.01 -28.0 -28.0 ... -11.1 -11.1 +Data variables: + elevation (latitude, longitude) float32 23MB -1.389e+03 ... -1.837e+03 +Attributes: + long_name: GDAL Band Number 1 + grid_mapping: crs
+
+
+
+
+
+
+
+In [ ]:
+
+
+
+
+
+dem.Dataset.to_netcdf("./data/iceland/dem.nc")
+
+
+
+
+
+
+
+
+
+Mesh¶
+
+
+
+
+
+
+
+In [ ]:
+
+
+
+
+
+mesh = pmesh.set(
+ mesh_generator='oceanmesh',
+ bgmesh = "om",
+ dem_source="./data/iceland/dem.nc",
+ type='tri2d',
+ geometry=geometry,
+ coastlines=coastlines,
+ grad = 0.15,
+ bathy_gradient= True,
+ resolution_min=0.02,
+ resolution_max=1.50,
+ alpha_wavelength= 100, # number of element to resolve WL
+ alpha_slope= 10, # number of element to resolve bathy gradient
+)
+mesh.Dataset.dims
+
+
+
+
+
+
+
+
+
+
+
+2024-07-15 07:22:44,770 WARNING pyposeidon.boundary __init__:67 coastlines is not a file, trying with geopandas Dataset +2024-07-15 07:22:46,870 INFO pyposeidon.moceanmesh get:157 Creating mesh with Oceanmesh + +2024-07-15 07:22:46,871 INFO pyposeidon.moceanmesh make_oceanmesh:420 Executing oceanmesh +2024-07-15 07:22:50,713 INFO pyposeidon.moceanmesh make_oceanmesh:464 oceanmesh: distance function +2024-07-15 07:22:50,978 INFO pyposeidon.moceanmesh make_oceanmesh:466 oceanmesh: feature size function +2024-07-15 07:22:53,083 INFO pyposeidon.moceanmesh make_oceanmesh:477 ./data/iceland/dem.nc +2024-07-15 07:22:53,191 INFO pyposeidon.moceanmesh make_oceanmesh:479 oceanmesh: WL size function +2024-07-15 07:22:53,556 INFO pyposeidon.moceanmesh make_oceanmesh:486 oceanmesh: bathymetric size function +2024-07-15 07:22:54,148 INFO pyposeidon.moceanmesh make_oceanmesh:539 oceanmesh: generate mesh +2024-07-15 07:23:16,682 INFO pyposeidon.moceanmesh make_oceanmesh:566 oceanmesh: boundaries ++
+
+Out[Â ]:
+
+
+FrozenMappingWarningOnValuesAccess({'nSCHISM_hgrid_node': 36683, 'nSCHISM_hgrid_face': 70876, 'nMaxSCHISM_hgrid_face_nodes': 3, 'bnodes': 2510})
+
+
+
+
+
+
+
+In [ ]:
+
+
+
+
+
+mesh.Dataset.pplot.mesh()
+
+
+
+
+
+
+
+
+Out[Â ]:
+
+
+<Axes: title={'center': 'Mesh plot'}, xlabel='Longitude (degrees)', ylabel='Latitude (degrees)'>
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+In [ ]:
+
+
+
+
+
+mesh.to_file('./data/iceland/hgrid.gr3')
+!head ./data/iceland/hgrid.gr3
+
+
+
+
+
+
+
+
+
+
+
+2024-07-15 07:23:17,325 INFO pyposeidon.mesh to_file:376 writing mesh to file ./data/iceland/hgrid.gr3 + uniform.gr3 + 70876 36683 +1 -21.5669000313 67.5843322685 0.0000000000 +2 -13.5993326425 64.9130533526 0.0000000000 +3 -23.9949065590 65.5776208081 0.0000000000 +4 -14.6220620374 64.3555872484 0.0000000000 +5 -17.3052891209 63.1876102642 0.0000000000 +6 -20.2997668446 63.4227245966 0.0000000000 +7 -18.2148689814 66.9584673197 0.0000000000 +8 -25.8802156161 63.4917894940 0.0000000000 ++
+
+
+
+
+
+
+
+
+Fix bathy¶
+
+
+
+
+
+
+
+In [ ]:
+
+
+
+
+
+#define in a dictionary the properties of the model..
+model_parameters = {
+ "solver_name": "telemac",
+ "rpath": "./data/iceland/",
+ "dem_source": "./data/iceland/dem.nc",
+ "mesh_file": "./data/iceland/hgrid.gr3",
+ "update": ["dem"], #set which component should be updated (meteo,dem,model)
+ "start_date": "2018-10-01T00:00:00",
+ "time_frame": "1d",
+ "global": False,
+}
+model_parameters
+
+
+
+
+
+
+
+
+Out[Â ]:
+
+
+{'solver_name': 'telemac',
+ 'rpath': './data/iceland/',
+ 'dem_source': './data/iceland/dem.nc',
+ 'mesh_file': './data/iceland/hgrid.gr3',
+ 'update': ['dem'],
+ 'start_date': '2018-10-01T00:00:00',
+ 'time_frame': '1d',
+ 'global': False}
+
+
+
+
+
+
+
+In [ ]:
+
+
+
+
+
+model = pmodel.set(**model_parameters)
+model.create()
+model.mesh.to_file('./data/iceland/hgrid.gr3')
+!head ./data/iceland/hgrid.gr3
+
+
+
+
+
+
+
+
+
+
+
+2024-07-15 07:23:18,386 INFO pyposeidon.tools open_dataset:350 extracting dem from ./data/iceland/dem.nc + +2024-07-15 07:23:18,387 DEBUG pyposeidon.tools open_dataset:361 Engine: netcdf4 +2024-07-15 07:23:18,393 WARNING pyposeidon.dem dem_:163 Lon must be within -28.006250000209224 and -11.097916666541362 +2024-07-15 07:23:18,393 WARNING pyposeidon.dem dem_:164 compensating if global dataset available +2024-07-15 07:23:18,394 WARNING pyposeidon.dem dem_:167 Lat is within 61.99374999984879 and 68.00624999975389 +2024-07-15 07:23:18,397 INFO pyposeidon.dem dem_:210 dem done + +2024-07-15 07:23:18,398 WARNING pyposeidon.dem __init__:81 coastlines not present, aborting adjusting dem + +2024-07-15 07:23:18,399 INFO pyposeidon.mesh read_file:239 read mesh file ./data/iceland/hgrid.gr3 +2024-07-15 07:23:18,734 INFO pyposeidon.telemac create:1007 bathy is set to 0, interpolating dem on mesh +2024-07-15 07:23:18,735 INFO pyposeidon.dem dem_on_mesh:224 .. interpolating on mesh .. + +2024-07-15 07:23:18,736 INFO pyposeidon.dem dem_on_mesh:233 .. using elevation values .. + ++
+
+
+
+
+/home/tomsail/miniconda3/envs/pos_test/lib/python3.11/site-packages/pyposeidon/mesh.py:247: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]` + ni, nj = df.iloc[0].str.split()[0] ++
+
+
+
+
+2024-07-15 07:23:19,344 INFO pyposeidon.dem dem_on_mesh:264 dem on mesh done + +2024-07-15 07:23:19,345 WARNING pyposeidon.telemac bath:871 coastlines not present, aborting adjusting dem + +2024-07-15 07:23:19,345 INFO pyposeidon.telemac get_depth_from_dem:887 Dem already adjusted + +2024-07-15 07:23:19,346 INFO pyposeidon.telemac get_depth_from_dem:898 updating bathymetry .. + +2024-07-15 07:23:19,347 INFO pyposeidon.telemac create:1009 found bathy in mesh file + +2024-07-15 07:23:19,348 INFO pyposeidon.telemac force:763 skipping meteo .. + +2024-07-15 07:23:19,349 INFO pyposeidon.telemac config:659 using default telemac2d json config file ... + +2024-07-15 07:23:19,345 WARNING pyposeidon.telemac bath:871 coastlines not present, aborting adjusting dem + +2024-07-15 07:23:19,345 INFO pyposeidon.telemac get_depth_from_dem:887 Dem already adjusted + +2024-07-15 07:23:19,346 INFO pyposeidon.telemac get_depth_from_dem:898 updating bathymetry .. + +2024-07-15 07:23:19,347 INFO pyposeidon.telemac create:1009 found bathy in mesh file + +2024-07-15 07:23:19,348 INFO pyposeidon.telemac force:763 skipping meteo .. + +2024-07-15 07:23:19,349 INFO pyposeidon.telemac config:659 using default telemac2d json config file ... + +2024-07-15 07:23:19,360 INFO pyposeidon.mesh to_file:376 writing mesh to file ./data/iceland/hgrid.gr3 + uniform.gr3 + 70876 36683 +1 -21.5669000313 67.5843322685 684.0000000000 +2 -13.5993326425 64.9130533526 66.0000000000 +3 -23.9949065590 65.5776208081 54.0000000000 +4 -14.6220620374 64.3555872484 23.0000000000 +5 -17.3052891209 63.1876102642 897.0000000000 +6 -20.2997668446 63.4227245966 10.0000000000 +7 -18.2148689814 66.9584673197 394.0000000000 +8 -25.8802156161 63.4917894940 355.0000000000 ++
+
+
+
+
+
+
+In [ ]:
+
+
+
+
+
+def parse_hgrid_nodes(path: os.PathLike[str] | str) -> pd.DataFrame:
+ with open(path, "rb") as fd:
+ _ = fd.readline()
+ _, no_points = map(int, fd.readline().strip().split(b" "))
+ content = io.BytesIO(b''.join(next(fd) for _ in range(no_points)))
+ nodes = pd.read_csv(
+ content,
+ engine="pyarrow",
+ sep="\t",
+ header=None,
+ names=["lon", "lat", "depth"],
+ index_col=0
+ )
+ nodes = nodes.reset_index(drop=True)
+ return nodes
+
+def parse_hgrid_elements3(path: os.PathLike[str] | str) -> pd.DataFrame:
+ with open(path, "rb") as fd:
+ _ = fd.readline()
+ no_elements, no_points = map(int, fd.readline().strip().split(b" "))
+ for _ in range(no_points):
+ next(fd)
+ content = io.BytesIO(b''.join(next(fd) for _ in range(no_elements)))
+ elements = pd.read_csv(
+ content,
+ engine="pyarrow",
+ sep="\t",
+ header=None,
+ names=["no_nodes", "n1", "n2", "n3"],
+ index_col=0
+ )
+ elements = elements.assign(
+ n1=elements.n1 - 1,
+ n2=elements.n2 - 1,
+ n3=elements.n3 - 1,
+ ).reset_index(drop=True)
+ return elements
+
+def get_skews_and_base_cfls(lons, lats, depths) -> np.ndarray:
+ # The shape of each one of the input arrays needs to be (3, <no_triangles>)
+ #ell = pymap3d.Ellipsoid.from_name("wgs84")
+ ell = pymap3d.Ellipsoid(6378206.4, 6378206.4, "schism", "schism")
+ local_x, local_y, _ = pymap3d.geodetic2enu(lats, lons, depths, lats[0], lons[0], depths[0], ell=ell)
+ areas = (local_x[1] * local_y[2] - local_x[2] * local_y[1]) * 0.5
+ rhos = np.sqrt(areas / np.pi)
+ max_sides = np.maximum(
+ np.sqrt(local_x[1] ** 2 + local_y[1] ** 2),
+ np.sqrt(local_x[2] ** 2 + local_y[2] ** 2),
+ np.sqrt((local_x[2] - local_x[1]) ** 2 + (local_y[2] - local_y[1]) ** 2),
+ )
+ skews = max_sides / rhos
+ base_cfls = np.sqrt(9.81 * np.maximum(0.1, depths.mean(axis=0))) / rhos / 2
+ return skews, base_cfls
+
+def get_skews_and_base_cfls_from_path(path: os.PathLike[str] | str) -> np.ndarray:
+ nodes = parse_hgrid_nodes(path)
+ elements = parse_hgrid_elements3(path)
+ tri = elements[["n1", "n2", "n3"]].values
+ lons = nodes.lon.values[tri].T
+ lats = nodes.lat.values[tri].T
+ depths = nodes.depth.values[tri].T
+ skews, base_cfls = get_skews_and_base_cfls(lons=lons, lats=lats, depths=depths)
+ return skews, base_cfls
+
+
+
+
+
+
+
+
+In [ ]:
+
+
+
+
+
+skews, base_cfls = get_skews_and_base_cfls_from_path("./data/iceland/hgrid.gr3")
+CFL_THRESHOLD = 0.4
+print(f"elements violating CFL threshold < {CFL_THRESHOLD}")
+print("time N %")
+for dt in (1, 50, 75, 100, 120, 150, 200, 300, 400, 600, 900, 1200, 1800, 3600):
+ violations = (base_cfls * dt < CFL_THRESHOLD).sum()
+ print(f"{dt:>4d} {violations:>12d} {violations / len(base_cfls) * 100:>8.2f}%")
+
+
+
+
+
+
+
+
+
+
+
+
+elements violating CFL threshold < 0.4 +time N % + 1 70876 100.00% + 50 4028 5.68% + 75 2096 2.96% + 100 1580 2.23% + 120 1377 1.94% + 150 1210 1.71% + 200 1029 1.45% + 300 825 1.16% + 400 777 1.10% + 600 48 0.07% + 900 0 0.00% +1200 0 0.00% +1800 0 0.00% +3600 0 0.00% ++
+
+
+
+
+
+
+In [ ]:
+
+
+
+
+
+pd.DataFrame({"skew": skews}).describe([0.25, 0.5, 0.75, 0.9, 0.95, 0.99, 0.995, 0.999])
+
+
+
+
+
+
+
+
+Out[Â ]:
+
+
+
+
+
+
+
+
+| + | skew | +
|---|---|
| count | +70876.000000 | +
| mean | +3.766559 | +
| std | +0.538883 | +
| min | +2.508742 | +
| 25% | +3.425276 | +
| 50% | +3.837101 | +
| 75% | +4.131190 | +
| 90% | +4.394984 | +
| 95% | +4.572409 | +
| 99% | +4.908128 | +
| 99.5% | +5.029456 | +
| 99.9% | +5.334526 | +
| max | +9.411692 | +
+
+
+
+
+
+
+In [ ]:
+
+
+
+
+
+def get_meta() -> gpd.GeoDataFrame:
+ meta_web = searvey.get_ioc_stations().drop(columns=["lon", "lat"])
+ meta_api = (
+ pd.read_json(
+ "http://www.ioc-sealevelmonitoring.org/service.php?query=stationlist&showall=all"
+ )
+ .drop_duplicates()
+ .drop(columns=["lon", "lat"])
+ .rename(columns={"Code": "ioc_code", "Lon": "longitude", "Lat": "latitude"})
+ )
+ merged = pd.merge(
+ meta_web,
+ meta_api[["ioc_code", "longitude", "latitude"]].drop_duplicates(),
+ on=["ioc_code"],
+ )
+ updated = merged.assign(
+ geometry=gpd.points_from_xy(merged.longitude, merged.latitude, crs="EPSG:4326")
+ )
+ return updated
+
+ioc_ = get_meta()
+ioc_[bbox.contains(ioc_.geometry)].to_csv('data/iceland/stations.csv')
+
+
+
+
+
+
+
+
+In [ ]:
+
+
+
+
+
+#define in a dictionary the properties of the model..
+model_parameters = {
+ "solver_name": "telemac",
+ "tag": "telemac2d",
+ "rpath": "./data/iceland/20181001",
+ "mesh_file": "./data/iceland/hgrid.gr3",
+ "update": ["all"], #set which component should be updated (meteo,dem,model)
+ "meteo_source": glob.glob("data/uvp_*.grib"),
+ "meteo_merge": "last", # combine meteo
+ "meteo_combine_by": "nested",
+ "meteo_xr_kwargs": {"concat_dim": "step"},
+ "start_date": "2018-10-01T00:00:00",
+ "time_frame": "2d",
+ "obs": "data/iceland/stations.csv",
+ "monitor": True,
+ "parameters": {
+ "dt": 100
+ }
+}
+model_parameters
+
+
+
+
+
+
+
+
+Out[Â ]:
+
+
+{'solver_name': 'telemac',
+ 'tag': 'telemac2d',
+ 'rpath': './data/iceland/20181001',
+ 'mesh_file': './data/iceland/hgrid.gr3',
+ 'update': ['all'],
+ 'meteo_source': ['data/uvp_2018100200.grib',
+ 'data/uvp_2018100112.grib',
+ 'data/uvp_2018100100.grib',
+ 'data/uvp_2018100212.grib'],
+ 'meteo_merge': 'last',
+ 'meteo_combine_by': 'nested',
+ 'meteo_xr_kwargs': {'concat_dim': 'step'},
+ 'start_date': '2018-10-01T00:00:00',
+ 'time_frame': '2d',
+ 'obs': 'data/iceland/stations.csv',
+ 'monitor': True,
+ 'parameters': {'dt': 100}}
+
+
+
+
+
+
+
+
+
+run days 1&2¶
+
+
+
+
+
+
+
+
+In [ ]:
+
+
+
+
+
+a = pmodel.set(**model_parameters)
+a.create()
+# IMPORTANT! Here, for a simple surge application,
+# we will need close all boundaries, otherwise the
+# model will run out of water
+a.mesh.Dataset.type[:] = 'closed' # it will create a cli file with all boundaries closed (this can be done only once)
+# we also need to drop some meteo variables, it is necesarry for zarr export
+a.output(**{"global": False})
+a.set_obs()
+a.save() # saves the json model file
+a.run() # runs the model
+
+
+
+
+
+
+
+
+
+
+
+2024-07-15 07:23:26,793 INFO pyposeidon.telemac create:987 no dem available + +2024-07-15 07:23:26,794 INFO pyposeidon.mesh read_file:239 read mesh file ./data/iceland/hgrid.gr3 ++
+
+
+
+
+/home/tomsail/miniconda3/envs/pos_test/lib/python3.11/site-packages/pyposeidon/mesh.py:247: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]` + ni, nj = df.iloc[0].str.split()[0] ++
+
+
+
+
+2024-07-15 07:23:27,266 INFO pyposeidon.telemac create:1009 found bathy in mesh file + +2024-07-15 07:23:27,267 INFO pyposeidon.meteo cfgrib:147 extracting meteo +2024-07-15 07:23:27,267 INFO pyposeidon.meteo cfgrib:147 extracting meteo +2024-07-15 07:23:32,795 INFO pyposeidon.meteo cfgrib:170 combining meteo by keeping the newest value +2024-07-15 07:23:32,818 INFO pyposeidon.meteo cfgrib:292 meteo done + +2024-07-15 07:23:32,820 INFO pyposeidon.telemac config:659 using default telemac2d json config file ... + +2024-07-15 07:23:32,822 INFO pyposeidon.telemac get_depth_from_dem:887 Dem already adjusted + +2024-07-15 07:23:32,822 INFO pyposeidon.telemac output:1039 Keeping bathymetry from hgrid.gr3 .. + +2024-07-15 07:23:32,824 INFO pyposeidon.telemac output:1041 saving geometry file.. +2024-07-15 07:23:32,831 DEBUG pyposeidon.telemac flip_triangles:170 reversed 0 triangles +2024-07-15 07:23:32,832 DEBUG pyposeidon.telemac flip_triangles:171 [TEMPORARY FIX]: REMOVE THE POLE TRIANGLES +2024-07-15 07:23:32,835 INFO pyposeidon.telemac mesh_to_slf:928 Manning file created.. + +2024-07-15 07:23:33,040 INFO pyposeidon.mesh to_file:376 writing mesh to file ./data/iceland/20181001/hgrid.gr3 +2024-07-15 07:23:34,071 INFO pyposeidon.telemac output:1048 saving meteo file.. +2024-07-15 07:23:43,555 INFO pyposeidon.telemac output:1071 saving boundary file.. +Islands +2024-07-15 07:23:44,499 INFO pyposeidon.telemac output:1076 saving CAS file.. +2024-07-15 07:23:44,515 INFO pyposeidon.telemac output:1080 output done + +2024-07-15 07:23:44,516 INFO pyposeidon.telemac set_obs:1384 get stations from data/iceland/stations.csv + +2024-07-15 07:23:44,522 INFO pyposeidon.telemac set_obs:1407 set in-situ measurements locations + +2024-07-15 07:23:44,547 INFO pyposeidon.telemac set_obs:1449 write out stations.csv file + +2024-07-15 07:23:44,549 INFO pyposeidon.telemac set_obs:1452 write out stations.in file + +2024-07-15 07:23:44,554 INFO pyposeidon.telemac run:1088 executing model + + ~> Checking keyword/rubrique coherence ++
+
+
+
+
+
+
+
+
+run days 3&4 from hotstart¶
+
+
+
+
+
+
+
+In [ ]:
+
+
+
+
+
+# restart model
+prev_ = pd.Timestamp('2018-10-01')
+next_ = pd.Timestamp('2018-10-03')
+end_ = pd.Timestamp('2018-10-05')
+ppath = os.path.join('data/iceland', prev_.strftime("%Y%m%d"))
+npath = os.path.join('data/iceland', next_.strftime("%Y%m%d"))
+m = pyposeidon.model.read(os.path.join(ppath, "telemac2d_model.json"))
+meteo = pmeteo.Meteo(glob.glob("data/uvp_*.grib"),meteo_merge= "last", meteo_combine_by= "nested", meteo_xr_kwargs= {"concat_dim": "step"},)
+rs = cast.set(
+ solver_name="telemac",
+ model=m,
+ ppath=ppath, # old path
+ cpath=npath, # new path
+ meteo=meteo.Dataset.sel(time=slice(next_, end_)).compute(),
+ sdate=next_, # new start date
+ end_date=end_, # new end date
+ start=next_, # start
+ copy=True,
+)
+b = rs.run(execute=False)
+
+
+
+
+
+
+
+
+
+
+
+2024-07-15 07:24:50,247 INFO pyposeidon.meteo cfgrib:147 extracting meteo +2024-07-15 07:24:54,405 INFO pyposeidon.meteo cfgrib:170 combining meteo by keeping the newest value +2024-07-15 07:24:54,415 INFO pyposeidon.meteo cfgrib:292 meteo done + +2024-07-15 07:24:54,562 DEBUG pyposeidon.utils.cast run:560 Copy model files +2024-07-15 07:24:54,564 DEBUG pyposeidon.utils.cast copy_files:51 copied src -> dst: /home/tomsail/work/python/pyPoseidon/Tutorial/data/iceland/20181001/geo.slf -> /home/tomsail/work/python/pyPoseidon/Tutorial/data/iceland/20181003/geo.slf +2024-07-15 07:24:54,565 DEBUG pyposeidon.utils.cast copy_files:51 copied src -> dst: /home/tomsail/work/python/pyPoseidon/Tutorial/data/iceland/20181001/geo.cli -> /home/tomsail/work/python/pyPoseidon/Tutorial/data/iceland/20181003/geo.cli +2024-07-15 07:24:54,566 DEBUG pyposeidon.utils.cast copy_files:51 copied src -> dst: /home/tomsail/work/python/pyPoseidon/Tutorial/data/iceland/20181001/station.in -> /home/tomsail/work/python/pyPoseidon/Tutorial/data/iceland/20181003/station.in +2024-07-15 07:24:54,566 DEBUG pyposeidon.utils.cast run:566 .. done +2024-07-15 07:24:54,567 DEBUG pyposeidon.utils.cast run:569 copy restart file +2024-07-15 07:24:54,567 INFO pyposeidon.utils.cast copy_or_symlink:68 Copying: /home/tomsail/work/python/pyPoseidon/Tutorial/data/iceland/20181001/restart_2D.slf -> /home/tomsail/work/python/pyPoseidon/Tutorial/data/iceland/20181003/prev_2D.slf +2024-07-15 07:24:54,573 INFO pyposeidon.utils.cast run:577 process meteo + +2024-07-15 07:24:54,574 INFO pyposeidon.telemac to_force:770 saving meteo file.. ++
+
+
+
+
+INFO:xarray_selafin.Serafin:Writing the output file: "/home/tomsail/work/python/pyPoseidon/Tutorial/data/iceland/20181003/input_wind.slf" ++
+
+
+
+
+2024-07-15 07:24:56,038 INFO pyposeidon.telemac config:655 reading parameter file /home/tomsail/work/python/pyPoseidon/Tutorial/data/iceland/20181001/telemac2d_model.json + +2024-07-15 07:24:56,039 INFO pyposeidon.telemac config:742 output telemac2d CAS file ... + +2024-07-15 07:24:56,050 INFO pyposeidon.telemac set_obs:1384 get stations from data/iceland/stations.csv + +2024-07-15 07:24:56,055 INFO pyposeidon.telemac set_obs:1407 set in-situ measurements locations + +2024-07-15 07:24:56,055 INFO pyposeidon.mesh read_file:239 read mesh file ./data/iceland/hgrid.gr3 ++
+
+
+
+
+/home/tomsail/miniconda3/envs/pos_test/lib/python3.11/site-packages/pyposeidon/mesh.py:247: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]` + ni, nj = df.iloc[0].str.split()[0] ++
+
+
+
+
+2024-07-15 07:24:56,406 INFO pyposeidon.telemac set_obs:1449 write out stations.csv file + +2024-07-15 07:24:56,417 INFO pyposeidon.telemac set_obs:1452 write out stations.in file + +2024-07-15 07:24:56,417 INFO pyposeidon.telemac set_obs:1452 write out stations.in file + +2024-07-15 07:24:56,419 INFO pyposeidon.utils.cast run:601 done for date : 20181003.00 ++
+
+
+
+
+
+
+
+In [ ]:
+
+
+
+
+
+b.run()
+
+
+
+
+
+
+
+
+
+
+
+2024-07-15 07:24:56,427 INFO pyposeidon.telemac run:1088 executing model + ++
+
+
+
+
+
+
+
+
+run 4 days model - check¶
+
+
+
+
+
+
+
+In [ ]:
+
+
+
+
+
+#define in a dictionary the properties of the model..
+model_parameters = {
+ "solver_name": "telemac",
+ "tag": "telemac2d",
+ "rpath": "./data/iceland/20181001-04",
+ "mesh_file": "./data/iceland/hgrid.gr3",
+ "update": ["all"], #set which component should be updated (meteo,dem,model)
+ "meteo_source": glob.glob("data/uvp_*.grib"),
+ "meteo_merge": "last", # combine meteo
+ "meteo_combine_by": "nested",
+ "meteo_xr_kwargs": {"concat_dim": "step"},
+ "start_date": "2018-10-01T00:00:00",
+ "time_frame": "4d",
+ "global": False,
+ "obs": "data/iceland/stations.csv",
+ "monitor": True,
+ "parameters": {
+ "dt": 100
+ }
+}
+model_parameters
+
+
+
+
+
+
+
+
+Out[Â ]:
+
+
+{'solver_name': 'telemac',
+ 'tag': 'telemac2d',
+ 'rpath': './data/iceland/20181001-04',
+ 'mesh_file': './data/iceland/hgrid.gr3',
+ 'update': ['all'],
+ 'meteo_source': ['data/uvp_2018100200.grib',
+ 'data/uvp_2018100112.grib',
+ 'data/uvp_2018100100.grib',
+ 'data/uvp_2018100212.grib'],
+ 'meteo_merge': 'last',
+ 'meteo_combine_by': 'nested',
+ 'meteo_xr_kwargs': {'concat_dim': 'step'},
+ 'start_date': '2018-10-01T00:00:00',
+ 'time_frame': '4d',
+ 'global': False,
+ 'obs': 'data/iceland/stations.csv',
+ 'monitor': True,
+ 'parameters': {'dt': 100}}
+
+
+
+
+
+
+
+
+In [ ]:
+
+
+
+
+
+c = pmodel.set(**model_parameters)
+c.create()
+# IMPORTANT! Here, for a simple surge application,
+# we will need close all boundaries, otherwise the
+# model will run out of water
+c.mesh.Dataset.type[:] = 'closed' # it will create a cli file with all boundaries closed (this can be done only once)
+# we also need to drop some meteo variables, it is necesarry for zarr export
+c.meteo.Dataset = c.meteo.Dataset.compute()
+c.output()
+c.set_obs()
+c.save() # saves the json model file
+c.run() # runs the model
+
+
+
+
+
+
+
+
+
+
+
+2024-07-15 07:25:58,336 INFO pyposeidon.telemac create:987 no dem available + +2024-07-15 07:25:58,337 INFO pyposeidon.mesh read_file:239 read mesh file ./data/iceland/hgrid.gr3 ++
+
+
+
+
+/home/tomsail/miniconda3/envs/pos_test/lib/python3.11/site-packages/pyposeidon/mesh.py:247: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]` + ni, nj = df.iloc[0].str.split()[0] ++
+
+
+
+
+2024-07-15 07:25:58,688 INFO pyposeidon.telemac create:1009 found bathy in mesh file + +2024-07-15 07:25:58,689 INFO pyposeidon.meteo cfgrib:147 extracting meteo +2024-07-15 07:26:03,008 INFO pyposeidon.meteo cfgrib:170 combining meteo by keeping the newest value +2024-07-15 07:26:03,019 INFO pyposeidon.meteo cfgrib:292 meteo done + +2024-07-15 07:26:03,020 INFO pyposeidon.telemac config:659 using default telemac2d json config file ... + +2024-07-15 07:26:03,333 INFO pyposeidon.telemac get_depth_from_dem:887 Dem already adjusted + +2024-07-15 07:26:03,333 INFO pyposeidon.telemac output:1039 Keeping bathymetry from hgrid.gr3 .. + +2024-07-15 07:26:03,335 INFO pyposeidon.telemac output:1041 saving geometry file.. +2024-07-15 07:26:03,339 DEBUG pyposeidon.telemac flip_triangles:170 reversed 0 triangles +2024-07-15 07:26:03,340 DEBUG pyposeidon.telemac flip_triangles:171 [TEMPORARY FIX]: REMOVE THE POLE TRIANGLES +2024-07-15 07:26:03,342 INFO pyposeidon.telemac mesh_to_slf:928 Manning file created.. + ++
+
+
+
+
+INFO:xarray_selafin.Serafin:Writing the output file: "./data/iceland/20181001-04/geo.slf" ++
+
+
+
+
+2024-07-15 07:26:03,502 INFO pyposeidon.mesh to_file:376 writing mesh to file ./data/iceland/20181001-04/hgrid.gr3 +2024-07-15 07:26:04,362 INFO pyposeidon.telemac output:1048 saving meteo file.. ++
+
+
+
+
+INFO:xarray_selafin.Serafin:Writing the output file: "./data/iceland/20181001-04/input_wind.slf" ++
+
+
+
+
+2024-07-15 07:26:06,961 INFO pyposeidon.telemac output:1071 saving boundary file.. +Islands +2024-07-15 07:26:07,951 INFO pyposeidon.telemac output:1076 saving CAS file.. +2024-07-15 07:26:07,959 INFO pyposeidon.telemac output:1080 output done + +2024-07-15 07:26:07,960 INFO pyposeidon.telemac set_obs:1384 get stations from data/iceland/stations.csv + +2024-07-15 07:26:07,964 INFO pyposeidon.telemac set_obs:1407 set in-situ measurements locations + +2024-07-15 07:26:08,000 INFO pyposeidon.telemac set_obs:1449 write out stations.csv file + +2024-07-15 07:26:08,003 INFO pyposeidon.telemac set_obs:1452 write out stations.in file + +2024-07-15 07:26:08,008 INFO pyposeidon.telemac run:1088 executing model + ++
+
+
+
+
+
+
+
+
+compare results¶
+
+
+
+
+
+
+
+In [ ]:
+
+
+
+
+
+res_2days = xr.open_dataset("data/iceland/20181001-04/results_2D.slf")
+res_day1 = xr.open_dataset("data/iceland/20181001/results_2D.slf")
+res_day2 = xr.open_dataset("data/iceland/20181003/results_2D.slf")
+
+
+
+
+
+
+
+
+In [ ]:
+
+
+
+
+
+node = 15000
+p1 = res_day1.S.isel(node = node).hvplot(label = "Day 1&2")
+p2 = res_day2.S.isel(node = node).hvplot(label = "Day 3&4")
+p3 = res_2days.S.isel(node = node).hvplot(label = "4 Days", line_dash='dashed')
+p1 * p2 * p3
+
+
+
+
+
+
+
+
+
+
+
+Out[Â ]:
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+compare with observations¶
+
+
+
+
+
+
+
+In [ ]:
+
+
+
+
+
+res_2days = xr.open_dataset("data/iceland/20181001-04/results_1D.slf")
+res_day1 = xr.open_dataset("data/iceland/20181001/results_1D.slf")
+res_day2 = xr.open_dataset("data/iceland/20181003/results_1D.slf")
+res_day2
+
+
+
+
+
+
+
+
+Out[Â ]:
+
+
+
+
+<xarray.Dataset> Size: 35kB
+Dimensions: (time: 1729, node: 1)
+Coordinates:
+ x (node) float32 4B ...
+ y (node) float32 4B ...
+ * time (time) datetime64[ns] 14kB 2018-10-03 ... 2018-10-05
+Dimensions without coordinates: node
+Data variables:
+ U (time, node) float32 7kB ...
+ V (time, node) float32 7kB ...
+ S (time, node) float32 7kB ...
+Attributes:
+ title: HISTORY FILE
+ language: en
+ float_size: 4
+ endian: >
+ params: (1, 0, 0, 0, 0, 0, 0, 1, 0, 1)
+ ipobo: [1]
+ ikle2: [[1]]
+ variables: {'U': ('VELOCITY U', 'M/S'), 'V': ('VELOCITY V', 'M/S'), 'S'...
+ date_start: (2018, 10, 1, 0, 0, 0)
+
+
+
+
+
+
+In [ ]:
+
+
+
+
+
+stations = pd.read_csv("data/iceland/stations.csv")
+stations
+
+
+
+
+
+
+
+
+Out[Â ]:
+
+
+
+
+
+
+
+
+| + | Unnamed: 0 | +ioc_code | +gloss_id | +country | +location | +connection | +contacts | +dcp_id | +last_observation_level | +last_observation_time | +... | +observations_arrived_per_month | +observations_expected_per_month | +observations_ratio_per_month | +observations_ratio_per_day | +sample_interval | +average_delay_per_day | +transmit_interval | +geometry | +longitude | +latitude | +
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | +1241 | +reyk | +229 | +Iceland | +Reykjavik | +ftp | +Icelandic Coast Guard, Hydrographic Department... | +NaN | +1.46 | +2020-03-20 12:28 | +... | +NaN | +NaN | +0 | +0 | +1' | +NaN | +1h | +POINT (-21.93333 64.15) | +-21.93333 | +64.15 | +
1 rows × 26 columns
+
+
+
+
+
+
+
+In [ ]:
+
+
+
+
+
+from searvey import ioc
+data = ioc.get_ioc_station_data('reyk', endtime="2018-10-05", period=30)
+data.index = data['time']
+data = data.drop(columns=['time'])
+data
+
+
+
+
+
+
+
+
+
+
+
+INFO:searvey.ioc:reyk: Retrieving data from: http://www.ioc-sealevelmonitoring.org/bgraph.php?code=reyk&output=tab&period=30&endtime=2018-10-05T00:00:00 ++
+
+Out[Â ]:
+
+
+
+
+
+
+
+
+| + | prs | +ioc_code | +
|---|---|---|
| time | ++ | + |
| 2018-09-05 00:00:00 | +2.727 | +reyk | +
| 2018-09-05 00:01:00 | +2.738 | +reyk | +
| 2018-09-05 00:02:00 | +2.744 | +reyk | +
| 2018-09-05 00:03:00 | +2.749 | +reyk | +
| 2018-09-05 00:04:00 | +2.752 | +reyk | +
| ... | +... | +... | +
| 2018-10-04 23:56:00 | +2.166 | +reyk | +
| 2018-10-04 23:57:00 | +2.178 | +reyk | +
| 2018-10-04 23:58:00 | +2.196 | +reyk | +
| 2018-10-04 23:59:00 | +2.214 | +reyk | +
| 2018-10-05 00:00:00 | +2.229 | +reyk | +
43120 rows × 2 columns
+
+
+
+
+
+
+
+In [ ]:
+
+
+
+
+
+# detide
+from analysea.tide import detide
+surge = detide(data["prs"],lat = 64.15)
+surge
+
+
+
+
+
+
+
+
+
+
+
+/home/tomsail/miniconda3/envs/pos_test/lib/python3.11/site-packages/utide/harmonics.py:16: RuntimeWarning: invalid value encountered in cast + nshallow = np.ma.masked_invalid(const.nshallow).astype(int) +/home/tomsail/miniconda3/envs/pos_test/lib/python3.11/site-packages/utide/harmonics.py:17: RuntimeWarning: invalid value encountered in cast + ishallow = np.ma.masked_invalid(const.ishallow).astype(int) - 1 ++
+
+
+
+
+solve: matrix prep ... solution ... done. ++
+
+Out[Â ]:
+
+
+time +2018-09-05 00:00:00 -0.011617 +2018-09-05 00:01:00 -0.007003 +2018-09-05 00:02:00 -0.007361 +2018-09-05 00:03:00 -0.008691 +2018-09-05 00:04:00 -0.011993 + ... +2018-10-04 23:56:00 -0.050843 +2018-10-04 23:57:00 -0.047389 +2018-10-04 23:58:00 -0.037939 +2018-10-04 23:59:00 -0.028494 +2018-10-05 00:00:00 -0.022052 +Name: prs, Length: 43120, dtype: float64+
+
+
+
+
+
+
+
+In [ ]:
+
+
+
+
+
+p1 = res_day1.S.isel(node = 0).hvplot(label = "Day 1")
+p2 = res_day2.S.isel(node = 0).hvplot(label = "Day 2")
+p3 = res_2days.S.isel(node = 0).hvplot(label = "2 Days", line_dash='dashed')
+obs_ = surge.loc["2018-10-01":"2018-10-05"].hvplot(label = "Observations", color = 'k', line_dash='dotted')
+p1 * p2 * p3 * obs_
+
+
+
+
+
+
+
+
+
+
+
+Out[Â ]:
+
+
+
+
+
+
+