Generate station.in per mesh

[pmav99]

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:

$ 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.

[pmav99]

For now I just added the stations.in for v1.2 and v2.2

[tomsail]

stations.in file needs to reside in the respective model (telemac or schism) folder.

For now, serialize is only supporting the schism format (see issue in pyposeidon), I'll add the TELEMAC format support to fix this.