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: crsIn [ ]:
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[Â ]: