closest_point_station.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import xarray as xr\n",
    "import requests\n",
    "from scipy.spatial import KDTree\n",
    "import numpy as np\n",
    "\n",
    "# Earth's radius in meters\n",
    "EARTH_RADIUS_METERS = 6371000\n",
    "\n",
    "def haversine(lon1, lat1, lon2, lat2):\n",
    "    \"\"\"Calculate the great-circle distance between two points on the Earth using the Haversine formula.\"\"\"\n",
    "    # Convert decimal degrees to radians\n",
    "    lon1, lat1, lon2, lat2 = map(np.radians, [lon1, lat1, lon2, lat2])\n",
    "    \n",
    "    # Haversine formula\n",
    "    dlon = lon2 - lon1\n",
    "    dlat = lat2 - lat1\n",
    "    a = np.sin(dlat / 2)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon / 2)**2\n",
    "    c = 2 * np.arcsin(np.sqrt(a))\n",
    "    distance = EARTH_RADIUS_METERS * c\n",
    "    return distance\n",
    "\n",
    "# Step 1: Load the NetCDF file and extract node coordinates\n",
    "netcdf_file = '/scratch/07174/soelem/global_15-20_highres/fort.63.nc'\n",
    "ds = xr.open_dataset(netcdf_file)\n",
    "\n",
    "# Extract longitude (x) and latitude (y) coordinates of the nodes\n",
    "lon = ds['x'].values\n",
    "lat = ds['y'].values\n",
    "coordinates = np.column_stack((lon, lat))  # Stack them as (lon, lat) pairs\n",
    "\n",
    "# Step 2: NOAA API to get tide gauge stations metadata (coordinates and station_id)\n",
    "def get_noaa_stations():\n",
    "    \"\"\"Retrieve all NOAA tide gauge stations.\"\"\"\n",
    "    base_url = 'https://api.tidesandcurrents.noaa.gov/mdapi/prod/webapi/stations.json'\n",
    "    response = requests.get(base_url)\n",
    "    if response.status_code == 200:\n",
    "        stations = response.json()['stations']\n",
    "        station_list = []\n",
    "        for station in stations:\n",
    "            try:\n",
    "                station_id = station['id']\n",
    "                lat = float(station['lat'])\n",
    "                lon = float(station['lng'])\n",
    "                station_list.append({'station_id': station_id, 'lat': lat, 'lon': lon})\n",
    "            except KeyError:\n",
    "                continue\n",
    "        return station_list\n",
    "    else:\n",
    "        raise Exception(f\"Error retrieving NOAA station data, status code {response.status_code}\")\n",
    "\n",
    "# Step 3: Find the closest station-node pairs and store the first 10 sorted by distance\n",
    "def find_closest_station_node_pairs(coordinates, stations, top_n=10):\n",
    "    \"\"\"Find the closest station-node pairs and return the top N sorted by distance using KDTree.\"\"\"\n",
    "    kdtree = KDTree(coordinates)\n",
    "    distances = []\n",
    "\n",
    "    for station in stations:\n",
    "        station_coords = (station['lon'], station['lat'])\n",
    "        distance, node_index = kdtree.query(station_coords)\n",
    "\n",
    "        # Convert the distance from degrees to meters using the Haversine formula\n",
    "        node_lon, node_lat = coordinates[node_index]\n",
    "        distance_meters = haversine(station['lon'], station['lat'], node_lon, node_lat)\n",
    "\n",
    "        # Store station ID, node index, and distance in meters\n",
    "        distances.append({\n",
    "            'station_id': station['station_id'],\n",
    "            'node_index': node_index,\n",
    "            'distance': distance_meters\n",
    "        })\n",
    "\n",
    "    # Sort the distances list by the 'distance' key\n",
    "    distances_sorted = sorted(distances, key=lambda x: x['distance'])\n",
    "\n",
    "    # Return the top N closest pairs\n",
    "    return distances_sorted[:top_n]\n",
    "\n",
    "# Get NOAA tide gauge stations\n",
    "stations = get_noaa_stations()\n",
    "\n",
    "# Step 4: Find and print the top 10 closest station-node pairs\n",
    "top_10_pairs = find_closest_station_node_pairs(coordinates, stations, top_n=20)\n",
    "\n",
    "# Print the top 10 pairs\n",
    "print(f\"{'Station ID':<12} {'Node Index':<12} {'Distance (meters)':<20}\")\n",
    "for pair in top_10_pairs:\n",
    "    print(f\"{pair['station_id']:<12} {pair['node_index']:<12} {pair['distance']:<20.2f}\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "era5",
   "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.9.19"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}