explore_WebCritech.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import requests\n",
    "import pandas as pd\n",
    "import hvplot.pandas\n",
    "import json\n",
    "import io\n",
    "import holoviews as hv\n",
    "import geoviews as gv\n",
    "import numpy as np\n",
    "import geopandas as gp\n",
    "import searvey\n",
    "from dateutil.relativedelta import *\n",
    "\n",
    "hv.extension('bokeh')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the base URL for the API\n",
    "sl_api  = \"https://webcritech.jrc.ec.europa.eu/SeaLevelsDb/api/Group/\"\n",
    "tad_api = \"https://webcritech.jrc.ec.europa.eu/TAD_server/api/Groups/Get\"\n",
    "\n",
    "api_device = \"https://webcritech.jrc.ec.europa.eu/api/Device/\"\n",
    "\n",
    "# Fetch the list of providers\n",
    "sl_response = requests.get(sl_api)\n",
    "tad_response = requests.get(tad_api)\n",
    "sl_data = json.loads(sl_response.text)  # Parse the JSON response\n",
    "tad_data = json.loads(tad_response.text)  # Parse the JSON response"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data availability"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### SeaLevelDB stations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize a list to hold all station data\n",
    "SL_data = []\n",
    "\n",
    "# Loop through the list of providers and fetch their stations\n",
    "for provider in sl_data:\n",
    "    provider_name = provider['Name']\n",
    "    stations_url = sl_api + provider_name + \"/Devices\"\n",
    "    stations_response = requests.get(stations_url)\n",
    "    stations_data = json.loads(stations_response.text)  # Parse the JSON response\n",
    "    # Loop through each station and extract data\n",
    "    for device in stations_data:\n",
    "        station_data = {\n",
    "            'Provider': provider_name,\n",
    "            'Id': device['Id'],\n",
    "            'Name': device['Name'],\n",
    "            'lat': device['Lat'],\n",
    "            'lon': device['Lon'],\n",
    "            'LastAccessStatus': device['CurrentStatus']['LastAccessStatus'],\n",
    "            'LastDate': device['CurrentStatus']['LastDate'],\n",
    "            'LastDate': device['CurrentStatus']['LastDate'],\n",
    "            'State': device['CurrentStatus']['State'],\n",
    "            'SyncStatus': device['CurrentStatus']['SyncStatus'],\n",
    "            'FileType': device['FileType'],\n",
    "            'GroupId': device['GroupId'],\n",
    "            'MovAvgNp': device['MovAvgNp'],\n",
    "            'Notes': device.get('Notes'),  # Use .get() to handle missing keys\n",
    "            'Source': device['Source']\n",
    "            # Add other fields as necessary\n",
    "        }\n",
    "        SL_data.append(station_data)\n",
    "SL_df = pd.DataFrame(SL_data)\n",
    "SL_df\n",
    "SL_df.to_csv('SeaLevelDb_stations.csv')\n",
    "SL_df.groupby(['Provider','State']).count()[['Id']].hvplot.bar(\n",
    "    rot = 90,\n",
    "    logy = True,\n",
    "    ylim = [0.5, 2000], \n",
    "    grid = True\n",
    ").opts(\n",
    "    width = 1200,\n",
    "    height = 800,\n",
    "    title = \"SeaLevelDb data availability\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "TAD_data Stations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "TAD_data = []\n",
    "# Loop through the list of providers and fetch their stations\n",
    "for provider in tad_data:\n",
    "    provider_name = provider['Name']\n",
    "    stations_url = f\"https://webcritech.jrc.ec.europa.eu/TAD_server/api/Groups/GetGeoJSON?group={provider_name}&maxLatency=10000000\"\n",
    "    stations_response = requests.get(stations_url)\n",
    "    stations_data = json.loads(stations_response.text)  # Parse the JSON response\n",
    "    if isinstance(stations_data, dict): stations_data = [stations_data]\n",
    "    # Loop through each station and extract data\n",
    "    if len(stations_data)>0:\n",
    "        for device in stations_data[0]['features']:\n",
    "            properties = device.get('properties', {})\n",
    "            geometry = device.get('geometry', {})\n",
    "            coordinates = geometry.get('coordinates', [None, None])\n",
    "            latency = properties.get('Latency', {}).get('Literal')\n",
    "            last_data_date_str = properties.get('LastData', {}).get('Date')\n",
    "            last_data_date = pd.Timestamp(last_data_date_str)\n",
    "            \n",
    "            # Determine if the measurement is within the last week\n",
    "            state = \"active\" if last_data_date and (pd.Timestamp.now() - last_data_date) <= pd.Timedelta(days=7) else \"inactive\"\n",
    "\n",
    "            station_data = {\n",
    "                'Provider': properties.get('Provider'),\n",
    "                'Id': device.get('id'),\n",
    "                'Name': properties.get('Name'),\n",
    "                'lat': coordinates[1],\n",
    "                'lon': coordinates[0],\n",
    "                'LastAccessStatus': properties.get('LastData', {}).get('Date'),\n",
    "                'LastDate': last_data_date,\n",
    "                'SyncStatus': latency,\n",
    "                'State': state,\n",
    "                'FileType': None,  # Update this if there is a corresponding field in the new data\n",
    "                'GroupId': None,  # Update this if there is a corresponding field in the new data\n",
    "                'MovAvgNp': None,  # Update this if there is a corresponding field in the new data\n",
    "                'Notes': properties.get('Notes'),  # Use .get() to handle missing keys\n",
    "                'Source': None\n",
    "                # Add other fields as necessary\n",
    "            }\n",
    "            TAD_data.append(station_data)\n",
    "TAD_df = pd.DataFrame(TAD_data)\n",
    "TAD_df\n",
    "TAD_df.to_csv('TAD_stations.csv')\n",
    "TAD_df.groupby(['Provider','State']).count()[['Id']].hvplot.bar(\n",
    "    rot = 90,\n",
    "    logy = True,\n",
    "    ylim = [0.5, 2000], \n",
    "    grid = True\n",
    ").opts(\n",
    "    width = 1600,\n",
    "    height = 800,\n",
    "    title = \"SeaLevelDb data availability\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Reference stations already in use in Seareport: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_stofs2d_meta():\n",
    "    stofs2d = pd.read_csv(\n",
    "        \"https://polar.ncep.noaa.gov/stofs/data/stofs_2d_glo_elev_stat_v2_1_0\",\n",
    "        names=[\"coords\", \"name\"],\n",
    "        sep=\"!\",\n",
    "        header=None,\n",
    "        skiprows=1\n",
    "    )\n",
    "    stofs2d = stofs2d.assign(\n",
    "        lon=stofs2d.coords.str.split(\"\\t\", n=1).str[0].astype(float),\n",
    "        lat=stofs2d.coords.str.strip().str.rsplit(\"\\t\", n=1).str[1].astype(float),\n",
    "        stofs2d_name=stofs2d.name.str.strip(),\n",
    "    ).drop(columns=[\"coords\", \"name\"])\n",
    "    return stofs2d\n",
    "\n",
    "\n",
    "def get_ioc_meta() -> gp.GeoDataFrame:\n",
    "    meta_web = searvey.get_ioc_stations().drop(columns=[\"lon\", \"lat\"])\n",
    "    meta_api = (\n",
    "        pd.read_json(\n",
    "            \"http://www.ioc-sealevelmonitoring.org/service.php?query=stationlist&showall=all\"\n",
    "        )\n",
    "        .drop_duplicates()\n",
    "        .drop(columns=[\"lon\", \"lat\"])\n",
    "        .rename(columns={\"Code\": \"ioc_code\", \"Lon\": \"lon\", \"Lat\": \"lat\"})\n",
    "    )\n",
    "    merged = pd.merge(\n",
    "        meta_web,\n",
    "        meta_api[[\"ioc_code\", \"lon\", \"lat\"]].drop_duplicates(),\n",
    "        on=[\"ioc_code\"],\n",
    "    )\n",
    "    return merged\n",
    "\n",
    "def get_coops_meta() -> gp.GeoDataFrame: \n",
    "    coops = searvey.get_coops_stations()\n",
    "    coops['lon'] = coops['geometry'].x\n",
    "    coops['lat'] = coops['geometry'].y\n",
    "    coops = coops.drop(columns='geometry')\n",
    "    return coops\n",
    "\n",
    "def merge_ioc_and_stofs(ioc: pd.DataFrame, stofs2d: pd.DataFrame, coops = pd.DataFrame) -> pd.DataFrame:\n",
    "    stations = pd.concat((ioc, stofs2d), ignore_index=True)\n",
    "    stations = stations.assign(unique_id=stations.ioc_code.combine_first(stations.stofs2d_name))\n",
    "    stations = pd.concat((stations, coops),ignore_index=True)\n",
    "    stations = stations.assign(unique_id=stations.unique_id.combine_first(stations.nws_id))\n",
    "    return stations\n",
    "\n",
    "ioc = get_ioc_meta()\n",
    "stofs2d = get_stofs2d_meta()\n",
    "coops = get_coops_meta()\n",
    "stations = merge_ioc_and_stofs(ioc=ioc, stofs2d=stofs2d, coops = coops.drop(columns=\"removed\"))\n",
    "stations['is_sat'] = stations.unique_id.str.contains('SA')\n",
    "stations['source'] = \"stofs\"\n",
    "stations.loc[~stations.ioc_code.isna(), 'source'] = \"ioc\"\n",
    "stations.loc[~stations.nws_id.isna(), 'source'] = \"coops\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare all stations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot1 = stations.drop(columns='geometry')[~stations['is_sat']].hvplot.points(\n",
    "    x = 'lon',\n",
    "    y='lat',\n",
    "    c = 'source',\n",
    "    cmap = 'fire',\n",
    "    geo = True,\n",
    "    s = 20,\n",
    "    legend = False\n",
    ")\n",
    "plot2 = pd.concat([SL_df,TAD_df]).hvplot.points(\n",
    "    x = 'lon',\n",
    "    y='lat',\n",
    "    c = 'Provider',\n",
    "    line_color = 'k',\n",
    "    line_width = 0.25,\n",
    "    cmap = 'glasbey',\n",
    "    geo = True,\n",
    "    s=150,\n",
    "    tiles=\"OSM\",\n",
    ")\n",
    "\n",
    "(plot2 * plot1).opts(\n",
    "    width = 1400,\n",
    "    height = 1070\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compare now with the data availability"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SL_df['LastDate'] =  SL_df.apply(lambda x:pd.Timestamp(x['LastDate']).tz_localize(None), axis = 1)\n",
    "TAD_df['LastDate'] = TAD_df.apply(lambda x:pd.Timestamp(x['LastDate']).tz_localize(None), axis = 1)\n",
    "SL_df.loc[SL_df['LastDate']>pd.Timestamp.now(), 'LastDate']= pd.Timestamp.now()\n",
    "SL_df.loc[SL_df['LastDate']<pd.Timestamp(2010,1,1), 'LastDate']= pd.Timestamp(2010,1,1)\n",
    "SL_df.LastDate.hvplot.hist(bins=100).opts(title = \"SeaLevelDb\")\n",
    "TAD_df.LastDate.hvplot.hist().opts(title = \"TAD Server\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SL_df['period_inactive'] = SL_df.apply(lambda x:(pd.Timestamp.now() - pd.Timestamp(x['LastDate']).tz_localize(None)).total_seconds()/3600/24/365.25, axis = 1)\n",
    "TAD_df['period_inactive'] = TAD_df.apply(lambda x:(pd.Timestamp.now() - pd.Timestamp(x['LastDate']).tz_localize(None)).total_seconds()/3600/24/365.25, axis = 1)\n",
    "\n",
    "all_stations = pd.concat([SL_df,TAD_df])\n",
    "all_stations.period_inactive.hvplot.hist().opts(title = \"years inactive\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_stations.to_csv('TAD_SLDB_stations.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_3 = SL_df.hvplot.points(\n",
    "    x = 'lon',\n",
    "    y='lat',\n",
    "    c = 'period_inactive',\n",
    "    s = 50,\n",
    "    geo = True,\n",
    "    cmap = 'fire_r',\n",
    "    tiles=\"OSM\",\n",
    ").opts(height = 550, title = 'years inactive')\n",
    "\n",
    "((plot2 * plot1).opts(height=550, title = 'Provider') + plot_3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finnish (FMI) and Indonesian stations (BIG) seem to be exclusive on WebCritech Server. \n",
    "\n",
    "Although they stopped recording data since January 2020. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_4 = TAD_df.hvplot.points(\n",
    "    x = 'lon',\n",
    "    y='lat',\n",
    "    c = 'period_inactive',\n",
    "    line_color = \"k\",\n",
    "    line_width = 0.25,\n",
    "    cmap = 'fire_r',\n",
    "    s = 50,\n",
    "    geo = True,\n",
    "    tiles=\"OSM\",\n",
    ").opts(title = \"years inactive\",height=550)\n",
    "\n",
    "((plot2 * plot1).opts(height=550, title = 'Provider') + plot_4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Most of TAD stations seem to have stopped recording in January 2024\n",
    "\n",
    "India (INCOIS) station seem to be exclusive on WebCritech server"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get Station detailed metadata (NOT FINISHED)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "we'd need to get the First Date information, which can be retrieve by scrapping the website: at the URL `https://webcritech.jrc.ec.europa.eu/SeaLevelsDb/Device/1014` (the device number is the `Provider` column) \n",
    "which is contained in `html.hqs-js > body > div > div.container > section.row > div#dev-last-details.col-md-5 > dl.dl-horizontal > dd` of the source HTML. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from bs4 import BeautifulSoup\n",
    "# A function to scrape the \"First Date\" from the given URL\n",
    "def get_first_date(device_number):\n",
    "    print('.', end='')\n",
    "    url = f\"https://webcritech.jrc.ec.europa.eu/SeaLevelsDb/Device/{device_number}\"\n",
    "    response = requests.get(url)\n",
    "    if response.status_code == 200:\n",
    "        soup = BeautifulSoup(response.content, 'html.parser')\n",
    "        # Find all dl-horizontal elements, which contain dt and dd pairs\n",
    "        dl_elements = soup.find_all('dl', class_='dl-horizontal')\n",
    "        for dl in dl_elements:\n",
    "            # Find all dt elements within the current dl element\n",
    "            dt_elements = dl.find_all('dt')\n",
    "            # print(dt_elements)\n",
    "            for index, dt in enumerate(dt_elements):\n",
    "                # Check if the dt element's text is 'First date'\n",
    "                if dt.get_text().strip() == 'First date':\n",
    "                    # Get the corresponding dd element that follows the dt element\n",
    "                    first_date_dd = dl.find_all('dd')[index]\n",
    "                    # print(first_date_dd.get_text().strip())\n",
    "                    # 1/0\n",
    "                    return first_date_dd.get_text().strip()\n",
    "    return None\n",
    "\n",
    "# Apply the function to each row and store the results in a new column\n",
    "SL_df['FirstDate'] = SL_df['Id'].apply(get_first_date)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from searvey.multi import multithread\n",
    "\n",
    "# rate_limit \n",
    "def fetch_SeaLevelDb_station(tmin: pd.Timestamp, tmax: pd.Timestamp, id: int, NSTEPS = 10000000):\n",
    "    url = f\"https://webcritech.jrc.ec.europa.eu/SeaLevelsDb/api/Device/{id}/Data?tMin={tmin.strftime('%Y-%m-%d')}&tMax={tmax.strftime('%Y-%m-%d')}&nPts={NSTEPS}&field=level&mode=CSV\"\n",
    "    stations_response = requests.get(url)\n",
    "    if stations_response.ok:\n",
    "        data = json.loads(stations_response.text)  # Use `loads` instead of `load`\n",
    "        df = pd.DataFrame(data)\n",
    "        return df\n",
    "    else:\n",
    "        print(f\"Error: {stations_response.status_code}\")\n",
    "        return pd.DataFrame()\n",
    "\n",
    "\n",
    "def fetch_TAD_station(tmin: pd.Timestamp, tmax: pd.Timestamp, id: int, NSTEPS = 10000000): \n",
    "    url = f\"https://webcritech.jrc.ec.europa.eu/TAD_server/api/Data/Get/{id}?tMin={tmin.strftime('%Y-%m-%d')}&tMax={tmax.strftime('%Y-%m-%d')}&nRec={NSTEPS}&mode=CSV\"\n",
    "    stations_response = requests.get(url)\n",
    "    if stations_response.ok:\n",
    "        data = io.StringIO(stations_response.text)  # Create a text stream object\n",
    "        df = pd.read_csv(data)  # Read the DataFrame from the text stream\n",
    "        return df\n",
    "    else:\n",
    "        print(f\"Error: {stations_response.status_code}\")\n",
    "        return pd.DataFrame()\n",
    "\n",
    "\n",
    "def get_TAD_stations(df:pd.DataFrame, device_var: str):\n",
    "    for s, station in df.iterrows():\n",
    "        print(s, station)\n",
    "        1/0\n",
    "\n",
    "\n",
    "def get_SLDB_stations(df:pd.DataFrame, device_var: str):\n",
    "    for s, station in df.iloc[::-1].iterrows():\n",
    "        print(station)\n",
    "        func_kwargs = []\n",
    "        for tstep in pd.date_range(pd.Timestamp(2020, 1, 1), station.LastDate, freq='MS'):\n",
    "            tmin = tstep\n",
    "            tmax = tstep + relativedelta(months=+1)\n",
    "            func_kwargs.append(\n",
    "                dict(\n",
    "                    tmin=tmin,\n",
    "                    tmax=tmax,\n",
    "                    Device=station[device_var],\n",
    "                ),\n",
    "            )\n",
    "        results = multithread(\n",
    "            func=fetch_SeaLevelDb_station,\n",
    "            func_kwargs=func_kwargs,\n",
    "            n_workers=10,\n",
    "            print_exceptions=False,\n",
    "            disable_progress_bar=False,\n",
    "        )   \n",
    "        for result in results:\n",
    "            if result.result is not None:\n",
    "                df = result.result\n",
    "                print(df)\n",
    "        1/0\n",
    "\n",
    "get_SLDB_stations(SL_df, \"Id\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "Device = 94\n",
    "url = f\"https://webcritech.jrc.ec.europa.eu/TAD_server/api/Data/Get/{Device}?tMin=2024-08-01&tMax=2024-09-01&nRec={10000000}&mode=CSV\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stations_response = requests.get(url)\n",
    "data = io.StringIO(stations_response.text)  # Create a text stream object\n",
    "df = pd.read_csv(data)  # Read the DataFrame from the text stream\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.hvplot.line( y='Lev RAD (m)')"
   ]
  }
 ],
 "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}