Validation_Protocol.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Protocol for skill assessment\n",
    "\n",
    "The period of the skill assessment is the whole year 2023: from 2023-01-01 to 2023-11-30. \n",
    "\n",
    "It is the most energetic period of the year for cylones, shown is the next figure ![figure](https://private-user-images.githubusercontent.com/18373442/311166159-412c182f-9dd5-43d7-936e-00393c6854f2.png?jwt=eyJhbGciOiJIUzI1NiIsInR5cCI6IkpXVCJ9.eyJpc3MiOiJnaXRodWIuY29tIiwiYXVkIjoicmF3LmdpdGh1YnVzZXJjb250ZW50LmNvbSIsImtleSI6ImtleTUiLCJleHAiOjE3MTM0MjkxODgsIm5iZiI6MTcxMzQyODg4OCwicGF0aCI6Ii8xODM3MzQ0Mi8zMTExNjYxNTktNDEyYzE4MmYtOWRkNS00M2Q3LTkzNmUtMDAzOTNjNjg1NGYyLnBuZz9YLUFtei1BbGdvcml0aG09QVdTNC1ITUFDLVNIQTI1NiZYLUFtei1DcmVkZW50aWFsPUFLSUFWQ09EWUxTQTUzUFFLNFpBJTJGMjAyNDA0MTglMkZ1cy1lYXN0LTElMkZzMyUyRmF3czRfcmVxdWVzdCZYLUFtei1EYXRlPTIwMjQwNDE4VDA4MjgwOFomWC1BbXotRXhwaXJlcz0zMDAmWC1BbXotU2lnbmF0dXJlPTE1ZjM2ZjZhMmQxZDdjNjUyZTJlOWU2MGYxZjJjNGE3NzcxZmFmNzBlOGE0OTAyZDY1ZTI0NTM5ZGU2MDgzOWQmWC1BbXotU2lnbmVkSGVhZGVycz1ob3N0JmFjdG9yX2lkPTAma2V5X2lkPTAmcmVwb19pZD0wIn0.onbyTy-cJv_DgcedHqPRQhbauxBuV2inXMvttdvYgZw)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "it seems that the begenning of the year was marked with cyclone in the south Indian ocean. \n",
    "\n",
    "Whereas during the fall, the cylones were more frequent in the Pacific and north Atlantic oceans. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1 - get information cyclone / event  \n",
    "We will extract the cylone track from the [IBtracks](https://www.ncei.noaa.gov/products/international-best-track-archive) data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "! mkdir -p data\n",
    "! wget https://www.ncei.noaa.gov/data/international-best-track-archive-for-climate-stewardship-ibtracs/v04r00/access/netcdf/IBTrACS.ALL.v04r00.nc -O data/IBTrACS.ALL.v04r00.nc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load IBtracks data \n",
    "import xarray as xr \n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import tqdm\n",
    "import glob\n",
    "import os\n",
    "import json\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.cm as cm\n",
    "from analysea.tide import detide \n",
    "import geopandas as gpd\n",
    "from scipy.spatial import cKDTree\n",
    "\n",
    "cyclones = xr.open_dataset('data/IBTrACS.ALL.v04r00.nc')\n",
    "cyclones"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# select only the summer 2023 period\n",
    "indice_storm_2023 = []\n",
    "def is_in_summer_2023(times):\n",
    "    test1 = [pd.Timestamp(t.decode()) < pd.Timestamp('2023-11-30') for t in times]\n",
    "    test2 = [pd.Timestamp(t.decode()) > pd.Timestamp('2023-01-01') for t in times]\n",
    "    test = np.logical_and(test1, test2)\n",
    "    return np.any(test)\n",
    "\n",
    "for i_storm in tqdm.tqdm(cyclones.storm.values[13300:]): # no need to loop over 19th/20th century!! \n",
    "    times = cyclones.isel(storm=i_storm).iso_time.values\n",
    "    if is_in_summer_2023(times): \n",
    "        indice_storm_2023.append(i_storm)\n",
    "# takes ~ 20sec to run"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def extract_specs(istorm): \n",
    "    event = cyclones.isel(storm=int(istorm))\n",
    "    x, y, time = np.array(event.lon), np.array(event.lat), np.array(event.time)\n",
    "    mask = ~np.isnan(x)\n",
    "    x = x[mask]\n",
    "    y = y[mask]\n",
    "    time = time[mask]\n",
    "    wind_max = np.array(event.usa_wind)[mask]\n",
    "    r64 = np.array(event.usa_rmw)[mask]/(111 * np.cos(np.deg2rad(y)))\n",
    "    name = np.array(event.name)\n",
    "    basin = np.array(event.basin)[mask]\n",
    "    return x, y, time, wind_max, r64, name, basin\n",
    "\n",
    "def plot_cyclone_track(x, y, wind_max, r64, name, ax = None):\n",
    "    if ax is None:\n",
    "        fig, ax = plt.subplots(1,1, figsize=(16,8))\n",
    "    im = ax.scatter(x, y, c = wind_max, label=name, s = 1, cmap = 'jet', vmin = 0, vmax = 100)\n",
    "    for ix, xi in enumerate(x):\n",
    "        circ = plt.Circle((x[ix], y[ix]), \n",
    "                          radius=r64[ix], \n",
    "                          fill=False, \n",
    "                          hatch = '////', \n",
    "                          color = cm.jet(np.min([1, wind_max[ix]/100])), \n",
    "                          alpha = 0.3)\n",
    "        ax.add_patch(circ)\n",
    "    # ax.axis('equal')\n",
    "    return ax, im"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# check 2023 cyclone tracks\n",
    "fig, ax = plt.subplots(1,1, figsize=(16,8))\n",
    "for istorm in indice_storm_2023:\n",
    "    x, y, time, wind_max, r64, name, basin = extract_specs(istorm)\n",
    "    ax, im = plot_cyclone_track(x, y, wind_max, r64, name, ax)\n",
    "    if max(wind_max) > 100: \n",
    "        print(f\"index IBtracks #{istorm}, name: {name}, max wind {np.max(wind_max)}, start: {np.min(time)} end: {np.max(time)}\")\n",
    "plt.colorbar(im, ax = ax, orientation = 'horizontal', fraction = 0.05, aspect = 60)\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2 - compare the cyclone tracks with tide gauges locations\n",
    "### 2.1 - get clean tide gauges \n",
    "from `ioc_cleanup` clean tide gauges: https://github.com/seareport/ioc_cleanup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get clean tide gauges \n",
    "def files_in_folder(folder, ext = '.csv'):\n",
    "    list_files = []\n",
    "    for item in glob.glob(folder + '*'):\n",
    "        tmp = item.split(ext)[0]\n",
    "        root, name = os.path.split(tmp)\n",
    "        if item.endswith(ext):\n",
    "            if ext == \".parquet\":\n",
    "                list_files.append(name)\n",
    "            elif ext == \".csv\":\n",
    "                list_files.append(name)\n",
    "    return list_files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "CLEAN_FOLDER = \"/home/tomsail/work/python/seareport_org/ioc_cleanup/clean/\"\n",
    "list_TG_2023 = files_in_folder(CLEAN_FOLDER, ext = '.parquet')\n",
    "len(list_TG_2023)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 - identify tide gauges with SEASET"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## get SEASET stations -- latest version\n",
    "seaset_full = pd.read_csv(\"/home/tomsail/Documents/work/python/oceanmodelling/seaset/Notebooks/catalog_full_updated.csv\", index_col=0)\n",
    "# some ioc code got remove for same location ex prin/prin2\n",
    "\n",
    "def is_similar_station(ioc_code, station_list):\n",
    "    return any(station.startswith(ioc_code) for station in station_list)\n",
    "\n",
    "ioc_seaset = seaset_full.dropna(subset='ioc_code')\n",
    "subset_2023 = ioc_seaset[ioc_seaset['ioc_code'].apply(is_similar_station, station_list=list_TG_2023)]\n",
    "subset_2023"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3 - refine to stations which are in the vicity of the cyclone tracks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dist_max = 5 # in degrees\n",
    "\n",
    "def dist(lon1, lat1, lon2, lat2):\n",
    "    return np.sqrt((lat2 - lat1) ** 2 + (lon2 - lon1) **2)\n",
    "\n",
    "def is_close_station(station, lons, lats, dist_max = 5):\n",
    "    lons[lons > 180] = lons[lons > 180] - 360\n",
    "    lon, lat = station\n",
    "    return np.any(dist(lon, lat, lons, lats) < dist_max)\n",
    "\n",
    "def get_tracks(cyclone_data):\n",
    "    lons = []\n",
    "    lats = []\n",
    "    for storm in cyclone_data.storm.values:\n",
    "        lons.extend(cyclone_data.isel(storm=storm).lon.values)\n",
    "        lats.extend(cyclone_data.isel(storm=storm).lat.values)\n",
    "    return lons, lats\n",
    "\n",
    "def subset_from_cyclone(df, cyclone_data, dist_max = 5):\n",
    "    lons, lats = get_tracks(cyclone_data)\n",
    "    close_stations = df[df\n",
    "    .apply(\n",
    "        lambda row: is_close_station(\n",
    "            (row['longitude'], row['latitude']), \n",
    "            np.array(lons), \n",
    "            np.array(lats), \n",
    "            dist_max = dist_max\n",
    "        ), \n",
    "        axis=1\n",
    "    )]\n",
    "    return close_stations\n",
    "\n",
    "def subset_from_tracks(df, lons, lats, dist_max = 5):\n",
    "    close_stations = df[df\n",
    "    .apply(\n",
    "        lambda row: is_close_station(\n",
    "            (row['longitude'], row['latitude']), \n",
    "            np.array(lons), \n",
    "            np.array(lats),\n",
    "            dist_max = dist_max\n",
    "        ), \n",
    "        axis=1\n",
    "    )]\n",
    "    return close_stations\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# superpose with 2023 cyclone tracks\n",
    "fig, ax = plt.subplots(1,1, figsize=(16,8))\n",
    "for istorm in indice_storm_2023:\n",
    "    x, y, time, wind_max, r64, name, basin = extract_specs(istorm)\n",
    "    ax, im = plot_cyclone_track(x, y, wind_max, r64, name, ax)\n",
    "stations_impacted_all = subset_from_cyclone(subset_2023, cyclones.isel(storm=indice_storm_2023))\n",
    "stations_impacted_all.plot.scatter(x='longitude', y='latitude', ax=ax, s= 100, c = 'r', marker = \"*\", edgecolor = 'k')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.4 - split the analysis into basins\n",
    "the basins in IBTracks are the following ones "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "basins = {\n",
    "    \"East Pacific\":  'EP',\n",
    "    \"North Atlantic\":'NA',\n",
    "    \"North Indian\":  'NI',\n",
    "    \"South Indian\":  'SI',\n",
    "    \"South Pacific\": 'SP',\n",
    "    \"West Pacific\":  'WP',\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "out = dict()\n",
    "for basin in basins.keys():\n",
    "    zone = basins[basin]\n",
    "    for istorm in indice_storm_2023:\n",
    "        x, y, time, wind_max, r64, name, basins_ = extract_specs(istorm)\n",
    "        list_ = [b.decode() for b in basins_]\n",
    "        if zone in list_:\n",
    "            stations_impacted = subset_from_tracks(subset_2023, x, y)\n",
    "            stations_impacted.plot.scatter(x='longitude', y='latitude',ax=ax, s= 100, c = 'r', marker = \"*\", edgecolor = 'k')\n",
    "            if len(stations_impacted) > 0:\n",
    "                params = {\n",
    "                    'name': str(name)[2:-1], #byte litteral\n",
    "                    'start' : pd.Timestamp(time[0]).strftime('%Y-%m-%d %H:%M:%S'),\n",
    "                    'end' : pd.Timestamp(time[-1]).strftime('%Y-%m-%d %H:%M:%S'),\n",
    "                    'stations': [c_ for c_ in stations_impacted.ioc_code.values]\n",
    "                }\n",
    "                out[str(istorm)] = params\n",
    "                print(params)\n",
    "with open('stations_impacted_2023.json', 'w') as f:\n",
    "    json.dump(out, f, indent=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3 - detide the selected stations\n",
    "create a folder for detide stations first"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "! mkdir -p data/surge"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SURGE_FOLDER = \"./data/surge/\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for ii,ioc_code in enumerate(subset_2023.ioc_code):\n",
    "    lat = subset_2023.iloc[ii].latitude\n",
    "    if not os.path.exists(SURGE_FOLDER + f\"{ioc_code}.parquet\"):\n",
    "        df = pd.read_parquet(CLEAN_FOLDER + f\"{ioc_code}.parquet\")\n",
    "        surge = detide(df[df.columns[0]], lat=lat, resample_detide = True)\n",
    "        surge.to_frame().to_parquet(SURGE_FOLDER + f\"{ioc_code}.parquet\")\n",
    "# 1 hour without resampling \n",
    "# 30 sec "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4 - compare with model results "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "BASE = \"/home/tomsail/Documents/work/python/pyPoseidon/Tutorial/models/\"\n",
    "v0 = \"/home/tomsail/work/python/pyPoseidon/Tutorial/models/v0.0/telemac/results_2D.nc\"\n",
    "v0p2 = \"/home/tomsail/work/python/pyPoseidon/Tutorial/models/v0.2/telemac/results_2D.nc\"\n",
    "# v2p0 = \"/home/tomsail/work/models/results/global-v2/202307_2D_tri.zarr\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## load result file -- 2D\n",
    "model_tel_v0 = xr.open_dataset(v0)\n",
    "model_tel_v0p2 = xr.open_dataset(v0p2)\n",
    "# model_tel_v2 = xr.open_dataset(v2p0)\n",
    "model_tel_v0p2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load json\n",
    "import json\n",
    "\n",
    "with open(\"./data/storms/cyclones_2023.json\", \"r\") as f:\n",
    "    events = json.load(f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "global variables for the following plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_2D = [model_tel_v0, model_tel_v0p2 ]\n",
    "version2D = [\"v0.0\", \"v0.2.0\",\"v0.2.1\", \"v2.0\"]\n",
    "alpha = [1, 0.6, 0.4, 0.3,  0.3]\n",
    "colors = ['blue', 'red', 'purple', 'brown', 'green']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "useful functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def is_overlapping(tris, meshx):\n",
    "    PIR = 180\n",
    "    x1, x2, x3 = meshx[tris].T\n",
    "    return np.logical_or(abs(x2 - x1) > PIR, abs(x3 - x1) > PIR)\n",
    "\n",
    "def extract_max_elev(ds, times): \n",
    "    maxH = ds.elev.sel(time=times, method='nearest').max(dim='time')\n",
    "    return maxH\n",
    "\n",
    "def plot_max_elev(ax, ds, times): \n",
    "    m = is_overlapping(ds.face_nodes,ds.longitude)\n",
    "    max_elev = extract_max_elev(ds,times)\n",
    "    im = ax.tricontourf(\n",
    "        ds.longitude,\n",
    "        ds.latitude,\n",
    "        ds.face_nodes[~m],\n",
    "        max_elev.values, \n",
    "        levels = np.arange(-0, 0.5, 0.02), \n",
    "        extend = 'both')\n",
    "    return im\n",
    "\n",
    "def closest_n_points(nodes, N, meshXY, dist_max=np.inf):\n",
    "    mytree = cKDTree(meshXY)\n",
    "    d_, indice = mytree.query(nodes, range(1, N + 1))\n",
    "    indice[d_ > dist_max] = -1\n",
    "    mask = indice != -1\n",
    "    return indice[mask].T, d_[mask].T\n",
    "\n",
    "def extract_t_elev_1D(ds, seaset_id):\n",
    "    idx_ds = np.where(ds.seaset_id == seaset_id)[0]\n",
    "    if len(idx_ds) > 0:\n",
    "        elev_ = ds.isel(node=idx_ds[0]).elev.values\n",
    "        t_ = [pd.Timestamp(ti) for ti in ds.isel(node=idx_ds[0]).time.values]\n",
    "    else: \n",
    "        print(f\"station: {ioc_code}, seaset_id: {seaset_id} not found in model\")\n",
    "        t_ = None; elev_ = None\n",
    "    return pd.Series(elev_, index=t_)\n",
    "\n",
    "def extract_t_elev_2D(ds, x, y):\n",
    "    lons, lats = ds.longitude.values, ds.latitude.values\n",
    "    indx, dist_ = closest_n_points(np.array([x, y]).T, 1, np.array([lons,lats]).T)\n",
    "    ds_ = ds.isel(node=indx[0])\n",
    "    elev_ = ds_.elev.values\n",
    "    t_ = [pd.Timestamp(ti) for ti in ds_.time.values]\n",
    "    return pd.Series(elev_, index=t_), np.round(dist_, 2)\n",
    "\n",
    "def get_corr(df1: pd.DataFrame, df2: pd.Series): \n",
    "    ts1, ts2 = df1.align(df2, axis = 0)\n",
    "    ts1 = ts1.interpolate()\n",
    "    nan_mask1 = pd.isna(ts1)\n",
    "    nan_mask2 = pd.isna(ts2)\n",
    "    nan_mask = np.logical_or(nan_mask1.values.T[0], nan_mask2.values)\n",
    "    ts1 = ts1[~nan_mask]\n",
    "    ts2 = ts2[~nan_mask]\n",
    "    corr = ts1.corr(ts2)\n",
    "    return np.round(corr, 2), ts1, ts2\n",
    "\n",
    "def get_percentiles(ts1, ts2):\n",
    "    x = np.arange(0, 0.99, 0.001)\n",
    "    x = np.hstack([x, np.arange(0.99, 1, 0.0001)])\n",
    "    pc1 = np.zeros(len(x))\n",
    "    pc2 = np.zeros(len(x))\n",
    "    for it, thd in enumerate(x):\n",
    "        pc1[it] = ts1.quantile(thd)\n",
    "        pc2[it] = ts2.quantile(thd)\n",
    "    return pc1, pc2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We do here comparison for each cyclone.\n",
    "\n",
    "For Europe there is the special case of [Babet](https://en.wikipedia.org/wiki/Storm_Babet), which is not in the IBTracks database.\n",
    "\n",
    "Autumn 2023 was an active season for northern europe ([source](https://en.wikipedia.org/wiki/2023%E2%80%9324_European_windstorm_season)), with important storms like [Ciaran](https://en.wikipedia.org/wiki/Storm_Ciar%C3%A1n) or [Babet](https://en.wikipedia.org/wiki/Storm_Babet). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib widget"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for storm in events:\n",
    "    fig, ax = plt.subplots(1,1, figsize=(14,8))\n",
    "    x, y, time, wind_max, r64, name, basin = extract_specs(storm)\n",
    "    xmin, xmax = np.min(x), np.max(x)\n",
    "    ymin, ymax = np.min(y), np.max(y)\n",
    "    tmin = pd.Timestamp(events[storm]['start'])\n",
    "    tmax = pd.Timestamp(events[storm]['end'])\n",
    "    name = events[storm]['name']\n",
    "    im = plot_max_elev(ax, ds_2D[-1], pd.date_range(tmin, tmax, freq='1h'))\n",
    "    ax, im2 = plot_cyclone_track(x, y, wind_max, r64, name, ax = ax)\n",
    "    ax.axis('equal')\n",
    "    if name == \"BABET\": xmin, xmax, ymin, ymax = -10, 20, 40, 70 #\n",
    "    ax.set_xlim(xmin - 5 , xmax + 5 )\n",
    "    ax.set_ylim(ymin - 5 , ymax + 5 )\n",
    "    ax.set_title(f\"cyclone {name}, id#{storm} from {pd.Timestamp(tmin)} to {pd.Timestamp(tmax)}\")\n",
    "    plt.colorbar(im2, ax = ax, orientation = 'horizontal', pad = 0.07, fraction = 0.05, aspect = 60, label = 'max wind speed (m/s)')\n",
    "    plt.colorbar(im, ax = ax, pad = 0.02, fraction = 0.05, label = 'max elevation (m)')\n",
    "    stations_impacted = subset_2023[subset_2023.ioc_code.isin(events[storm]['stations'])]\n",
    "    plt.tight_layout()\n",
    "    # \n",
    "    fig1, ax1 = plt.subplots(len(stations_impacted), 1, figsize=(14,3*len(stations_impacted)))\n",
    "    if len(stations_impacted) == 1: ax1 = [ax1]\n",
    "    for i_s, ioc_code in enumerate(stations_impacted.ioc_code):\n",
    "        s = stations_impacted.iloc[i_s]\n",
    "        xl, yl = s.longitude, s.latitude\n",
    "        ax.scatter(xl, yl,\n",
    "            s=100, \n",
    "            lw=0.5, \n",
    "            c = colors[i_s % len(colors)], \n",
    "            marker = \"*\", \n",
    "            edgecolors = 'white', \n",
    "            label = ioc_code)\n",
    "        obs = pd.read_parquet(SURGE_FOLDER + f\"{ioc_code}.parquet\")\n",
    "        seaset_id = s.seaset_id\n",
    "        \n",
    "        # observations\n",
    "        obs = obs.loc[tmin:tmax]\n",
    "        ax1[i_s].plot(obs.index, obs.values , label=f\"{ioc_code}\", color = 'k', linestyle = '--')\n",
    "\n",
    "        # models\n",
    "        for ids, ds in enumerate(ds_2D):\n",
    "            mod , d_= extract_t_elev_2D(ds, xl, yl)\n",
    "            corr, ts1, ts2 = get_corr(mod, obs[obs.columns[0]])\n",
    "            ax1[i_s].plot(mod.index,mod.values , label=f\"model 2D {version2D[ids]}\", color = colors[i_s % len(colors)], alpha = alpha[ids])\n",
    "\n",
    "        ax1[i_s].set_xlim(tmin, tmax)\n",
    "        ax1[i_s].legend()\n",
    "    ax1[0].set_title(f\"cyclone {name}, id#{storm}\")\n",
    "    ax.legend()\n",
    "    plt.tight_layout()\n",
    "    print(storm, name)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "what about the for the whole period of simulation?\n",
    "\n",
    "Some storms might not have been saved in the IBTracks database but still induce a significant surge in some stations.\n",
    "\n",
    "Let's check the signal for all the stations in the surge folder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmin = pd.Timestamp(2023,1,1)\n",
    "tmax = pd.Timestamp(2023,10,31)\n",
    "\n",
    "list_surge_manual = ['viti', 'prus', 'mill', 'pkem', 'vard', 'wood', 'kahu', 'cald2', 'hmda', 'kinl', 'benc', 'dzao', \n",
    "                     'lerw2', 'live', 'cres', 'hie2', 'ptmt', 'ilfa', 'wick', 'rorv', 'pwil2', 'dkwa', 'bame', 'abed', 'lime', 'nshi', \n",
    "                     'newl',  'chst', 'bamf', 'kush', 'greg', 'bamd', 'honn', 'dutc',  'talc2', 'dpnc', 'quir2', 'fpnt', 'chrp', 'ishig', \n",
    "                     'bgct',  'vhbc', 'thev', 'npor', 'pslu', 'cher', 'boma', 'stjo', 'djve', 'ross', 'elak', 'hana', 'saig', 'pich2', \n",
    "                     'dapi', 'ptal2', 'stor', 'ande', 'aren', 'ohig3', 'newl2', 'trst', 'kusm', 'alam', 'leit', 'sado', 'asto', 'shee', \n",
    "                     'malo', 'corr2', 'prin2', 'nkfa', 'heys', 'helg', 'kawa', 'atka', 'qtro2', 'dove', 'waka', 'treg', 'cuxh', 'tosa', \n",
    "                     'yaku', 'brom', 'smog', 'darw', 'cwfl', 'huat', 'sprg', 'fue2', 'abur', 'guam', 'kwfl', 'pcha2', 'barn', \n",
    "                     'bapj', 'amal', 'hilo', 'wpwa', 'plym', 'coru', 'gokr', 'pmur', 'whit', 'nhav', 'porp', 'kungr', 'mumb', 'stqy2', 'crom', \n",
    "                     'sitk', 'stqy', 'oste', 'acnj', 'bang', 'naga', 'mare', 'fuka', 'herb', 'pagb', 'work', 'mhav', 'lajo', 'harw', \n",
    "                     'omae', 'coqu2', 'holy2', 'horn', 'sdpt', 'lowe', 'naha']\n",
    "surge_stations = subset_2023[subset_2023.ioc_code.isin(list_surge_manual)]\n",
    "countries = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))\n",
    "\n",
    "for i_s, ioc_code in tqdm.tqdm(enumerate(surge_stations.ioc_code)):        \n",
    "    # get station coordinates\n",
    "    s = surge_stations.iloc[i_s]\n",
    "    xl, yl = s.longitude, s.latitude\n",
    "    # get observations\n",
    "    obs = pd.read_parquet(SURGE_FOLDER + f\"{ioc_code}.parquet\")\n",
    "    obs = obs.loc[tmin:tmax]\n",
    "    seaset_id = s.seaset_id\n",
    "    # \n",
    "    mod , d_= extract_t_elev_2D(ds_2D[-1], xl, yl) # we take the v0.2 results only\n",
    "    corr, ts2, ts1 = get_corr(mod, obs[obs.columns[0]])\n",
    "    if corr > 0.0: \n",
    "        # if coefficient is too low, we don't plot it\n",
    "        fig = plt.figure(figsize=(16,5))\n",
    "        gs = fig.add_gridspec(1,3,  width_ratios=(2, 1, 1),\n",
    "                            left=0, right=0.99, bottom=0, top=0.99,\n",
    "                            wspace=0.07, hspace=0.03)\n",
    "        ax_plot1 = fig.add_subplot(gs[ 0])\n",
    "        ax_plot2 = fig.add_subplot(gs[ 1], sharey = ax_plot1)\n",
    "        ax_map = fig.add_subplot(gs[ 2])\n",
    "\n",
    "        ax.scatter(xl, yl,\n",
    "            s=100, \n",
    "            lw=0.5, \n",
    "            c = colors[i_s % len(colors)], \n",
    "            marker = \"*\", \n",
    "            edgecolors = 'white', \n",
    "            label = ioc_code)\n",
    "        \n",
    "        # models\n",
    "        for ids, ds in enumerate(ds_2D[-1:]):\n",
    "            mod , d_= extract_t_elev_2D(ds, xl, yl)\n",
    "            corr, ts2, ts1 = get_corr(mod, obs[obs.columns[0]])\n",
    "            mod95 = mod.quantile(0.95)\n",
    "            corr95, _, _ = get_corr(mod[mod > mod95], obs[obs.columns[0]])\n",
    "            mod99 = mod.quantile(0.99)\n",
    "            corr99, _, _ = get_corr(mod[mod > mod99], obs[obs.columns[0]])\n",
    "            ax_plot1.plot(mod.index,mod.values , label=f\"model {version2D[ids]}\", color = colors[i_s % len(colors)], alpha = alpha[ids])\n",
    "        ax_plot2.scatter(ts1.values, ts2.values, c= 'k', label=f\"model {version2D[ids]}, dist={d_}, Cr={corr}\", s=1, alpha = 0.3)\n",
    "        ax_plot2.set_xlabel('measured data')\n",
    "        ax_plot2.set_ylabel('modelled data')\n",
    "        pc1, pc2 = get_percentiles(ts1, ts2)\n",
    "        ax_plot2.axline([0,0],slope = 1, lw =1, linestyle = '--', color = 'k')\n",
    "        ax_plot2.scatter(pc1, pc2, color = colors[i_s % len(colors)], alpha = alpha[ids])\n",
    "        ax_plot2.plot(pc1, pc2, color = colors[i_s % len(colors)], alpha = alpha[ids])\n",
    "        # observations\n",
    "        ax_plot1.plot(obs.index, obs.values , label=f\"{ioc_code}\", color = 'k', linestyle = '--')\n",
    "        ax_plot1.set_xlim(tmin, tmax)\n",
    "        ax_plot1.legend()\n",
    "        ax_plot2.legend()\n",
    "        ax_plot1.grid(axis='both', color = 'grey')\n",
    "        ax_plot2.grid(axis='both', color = 'grey')\n",
    "        # map\n",
    "        ax_map.scatter(xl, yl, marker = \"*\", c = 'r')\n",
    "        _ = countries.plot(color='lightgrey', ax=ax_map, zorder=-1)\n",
    "        ax_map.set_xlim(xl-20, xl+20)\n",
    "        ax_map.set_ylim(yl-20, yl+20)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".venv",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}