Open pmav99 opened 1 month ago
For each mesh we should also create a corresponding station.in.
station.in
The following snippet is merging STOFS-2D and IOC stations in a single file. The results is similar to this:
$ head station.in 1 0 0 0 0 0 0 0 0 ! https://schism-dev.github.io/schism/master/input-output/optional-inputs.html#stationin-bp-format 3323 ! number of stations 1 144.2806072235 44.0247266972 0 ! abas 2269385 144.2900000000 44.0200000000 13.0000000000 916.6594922518 2 -2.0464618098 57.1249802456 0 ! abed 962967 -2.0800000000 57.1400000000 24.0000000000 2623.9968649481 3 131.4139218534 31.5887292407 0 ! abur 1854768 131.4100000000 31.5800000000 26.0000000000 1039.3086788890 4 -89.8377180390 13.5701359897 0 ! acaj 2746534 -89.8381280000 13.5737920000 9.0000000000 408.9377541903 5 -99.9228904646 16.8436445908 0 ! acap 2681197 -99.9166000000 16.8333000000 15.0000000000 1330.9068274144 6 -99.9228904646 16.8436445908 0 ! acap2 2681197 -99.9030000000 16.8379333000 15.0000000000 2210.0747947232 7 -74.4102481571 39.3729424866 0 ! acnj 1727855 -74.4183000000 39.3550000000 -0.0000000000 2111.7822884706 8 -74.4102481571 39.3729424866 0 ! acnj2 1727855 -74.4183000000 39.3550000000 -0.0000000000 2111.7822884706 $ tail station.in 3314 -133.2311735220 54.4803069628 0 ! "UJ845 SOUS00 SA845" 2077033 -133.2200000000 54.4600000000 345.0000000000 2370.6561082907 3315 -56.6849517955 54.4792339399 0 ! "UJ846 SOUS00 SA846" 2014221 -56.6800000000 54.4600000000 103.0000000000 2162.5208329908 3316 -36.8209103440 54.4617331441 0 ! "UJ847 SOUS00 SA847" 2527088 -36.8400000000 54.4600000000 2839.0000000000 1248.7824430343 3317 -17.1246075148 54.5407370460 0 ! "UJ848 SOUS00 SA848" 492125 -17.0000000000 54.4600000000 2442.0000000000 12055.4642648129 3318 147.4342946221 58.9668723570 0 ! "UJ849 SOUS00 SA849" 2908985 147.4100000000 58.9600000000 111.0000000000 1588.6814030428 3319 -178.5729671215 58.9347728367 0 ! "UJ850 SOUS00 SA850" 1232699 -178.5800000000 58.9600000000 2763.0000000000 2833.9881290880 3320 -144.4913277716 58.9552577946 0 ! "UJ851 SOUS00 SA851" 474508 -144.5700000000 58.9600000000 3772.0000000000 4541.7870224183 3321 -68.0504057985 58.9667668588 0 ! "UJ852 SOUS00 SA852" 1587142 -68.0200000000 58.9600000000 48.0000000000 1898.6426578148 3322 -50.7457083278 58.9051883153 0 ! "UJ853 SOUS00 SA853" 2699858 -51.0100000000 58.9600000000 3433.0000000000 16344.3736645492 3323 -25.5662962363 58.9674260336 0 ! "UJ854 SOUS00 SA854" 2101694 -25.5000000000 58.9600000000 2517.0000000000 3889.4252078064
import os import geopandas as gpd import holoviews as hv import hvplot.pandas import numpy as np import pandas as pd import panel as pn import searvey import sklearn.neighbors from pyposeidon.utils.cfl import parse_hgrid_nodes from pyposeidon.utils.cpoint import find_nearest_nodes from pyposeidon.utils.obs import serialize_stations EARTH_RADIUS = 6371000 # in meters @pn.cache def get_stofs2d_meta(): stofs2d = pd.read_csv( "https://polar.ncep.noaa.gov/stofs/data/stofs_2d_glo_elev_stat_v2_1_0", names=["coords", "name"], sep="!", header=None, skiprows=1 ) stofs2d = stofs2d.assign( lon=stofs2d.coords.str.split("\t", n=1).str[0].astype(float), lat=stofs2d.coords.str.strip().str.rsplit("\t", n=1).str[1].astype(float), stofs2d_name=stofs2d.name.str.strip(), ).drop(columns=["coords", "name"]) return stofs2d @pn.cache def get_ioc_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": "lon", "Lat": "lat"}) ) merged = pd.merge( meta_web, meta_api[["ioc_code", "lon", "lat"]].drop_duplicates(), on=["ioc_code"], ) return merged def find_closest_nodes(mesh_nodes: pd.DataFrame, points: pd.DataFrame): # The only requirement is that both `mesh_nodes and `points` have `lon/lat` columns tree = sklearn.neighbors.BallTree( np.radians(mesh_nodes[["lat", "lon"]]).values, metric="haversine", ) distances, indices = tree.query(np.radians(points[["lat", "lon"]].values)) closest_nodes = ( mesh_nodes .rename(columns={"lon": "mesh_lon", "lat": "mesh_lat"}) .iloc[indices.flatten()] .assign(distance=(distances.flatten() * EARTH_RADIUS)) .reset_index(names=["mesh_index"]) ) return pd.concat((points, closest_nodes), axis="columns") def merge_ioc_and_stofs(ioc: pd.DataFrame, stofs2d: pd.DataFrame) -> pd.DataFrame: stations = pd.concat((ioc, stofs2d), ignore_index=True) stations = stations.assign(unique_id=stations.ioc_code.combine_first(stations.stofs2d_name)) return stations mesh_nodes = parse_hgrid_nodes("/home/panos/Prog/poseidon/seareport_meshes/v2.2/schism/hgrid.gr3") ioc = get_ioc_meta() stofs2d = get_stofs2d_meta() stations = merge_ioc_and_stofs(ioc=ioc, stofs2d=stofs2d) closest_nodes = find_nearest_nodes(mesh_nodes=mesh_nodes, points=stations) serialize_stations(closest_nodes, path="./station.in")
Some of the functions live in the dev branch of pyposeidon but we could easily inline them.
dev
For now I just added the stations.in for v1.2 and v2.2
For each mesh we should also create a corresponding
station.in
.The following snippet is merging STOFS-2D and IOC stations in a single file. The results is similar to this:
Some of the functions live in the
dev
branch of pyposeidon but we could easily inline them.