RMSEv3.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "from scipy.spatial import cKDTree\n",
    "import netCDF4\n",
    "from noaa_coops import Station\n",
    "from sklearn.metrics import mean_squared_error\n",
    "from math import sqrt\n",
    "import time\n",
    "\n",
    "# Function to calculate RMSE\n",
    "def calculate_rmse(true_values, predicted_values):\n",
    "    if len(true_values) > 384 and len(predicted_values) > 384:\n",
    "        true_values = true_values[384:]\n",
    "        predicted_values = predicted_values[384:]\n",
    "    else:\n",
    "        return np.nan\n",
    "\n",
    "    valid_mask = ~np.isnan(true_values) & ~np.isnan(predicted_values)\n",
    "    true_values, predicted_values = true_values[valid_mask], predicted_values[valid_mask]\n",
    "\n",
    "    return sqrt(mean_squared_error(true_values, predicted_values))\n",
    "\n",
    "# Read the station details CSV file\n",
    "start_time = time.time()\n",
    "station_details = pd.read_csv('station_details.csv')\n",
    "print(f\"Reading station details took {time.time() - start_time:.2f} seconds\")\n",
    "\n",
    "# Load the NetCDF file\n",
    "start_time = time.time()\n",
    "nc_file = netCDF4.Dataset('/scratch/07174/soelem/test_sentinel_highres/fort.63.nc', 'r')\n",
    "print(f\"Loading NetCDF file took {time.time() - start_time:.2f} seconds\")\n",
    "\n",
    "latitudes = nc_file.variables['y'][:]\n",
    "longitudes = nc_file.variables['x'][:]\n",
    "\n",
    "# Normalize longitudes from 0-360 to -180 to +180 if needed\n",
    "longitudes = np.where(longitudes > 180, longitudes - 360, longitudes)\n",
    "\n",
    "# Create an array of tuples with the latitude and longitude\n",
    "points = np.column_stack((latitudes, longitudes))\n",
    "\n",
    "# Create a KDTree for quick nearest-neighbor lookup\n",
    "start_time = time.time()\n",
    "kdtree = cKDTree(points)\n",
    "print(f\"Building KDTree took {time.time() - start_time:.2f} seconds\")\n",
    "\n",
    "# Find the closest points in the NetCDF data for each station\n",
    "start_time = time.time()\n",
    "closest_points = []\n",
    "for index, row in station_details.iterrows():\n",
    "    _, closest_point_idx = kdtree.query([row['Latitude'], row['Longitude']])\n",
    "    closest_points.append(closest_point_idx)\n",
    "print(f\"Finding closest points took {time.time() - start_time:.2f} seconds\")\n",
    "\n",
    "# Print the array of closest points\n",
    "print(\"Closest points to stations:\", closest_points)\n",
    "\n",
    "# Add the closest point indices to the station details\n",
    "station_details['node'] = closest_points\n",
    "\n",
    "# Prepare a list to hold RMSE values\n",
    "rmse_values = []\n",
    "\n",
    "# Load the zeta values corresponding to the closest points only\n",
    "start_time = time.time()\n",
    "for index, row in station_details.iterrows():\n",
    "    try:\n",
    "        station_id = row['Station ID']\n",
    "        closest_point_idx = closest_points[index]\n",
    "\n",
    "        # Retrieve 'zeta' variable data for the closest point, skipping the first 384 values\n",
    "        nc_zeta = nc_file.variables['zeta'][384:, closest_point_idx]\n",
    "\n",
    "        # Initialize the station object\n",
    "        station = Station(station_id)\n",
    "\n",
    "        try:\n",
    "            station_data = station.get_data(\n",
    "                begin_date=\"20230716\",\n",
    "                end_date=\"20231031\",\n",
    "                product=\"hourly_height\",\n",
    "                datum=\"MSL\",\n",
    "                units=\"metric\",\n",
    "                time_zone=\"lst\"\n",
    "            )\n",
    "        except Exception as e:\n",
    "            print(f\"Error retrieving data for station {station_id}: {e}\")\n",
    "            rmse_values.append(np.nan)\n",
    "            continue\n",
    "\n",
    "        station_df = pd.DataFrame(station_data)\n",
    "\n",
    "        # Adjusted starting date after skipping 384 hours (16 days)\n",
    "        nc_zeta_series = pd.Series(nc_zeta.flatten(), index=pd.date_range(start=\"2023-08-01 00:00:00\", periods=len(nc_zeta), freq='H'))\n",
    "        station_df_resampled = station_df['v'].reindex(nc_zeta_series.index, method='nearest')\n",
    "\n",
    "        rmse = calculate_rmse(station_df_resampled.values, nc_zeta_series.values)\n",
    "        rmse_values.append(rmse)\n",
    "\n",
    "    except Exception as e:\n",
    "        print(f\"An unexpected error occurred for station {station_id}: {e}\")\n",
    "        rmse_values.append(np.nan)\n",
    "        continue\n",
    "\n",
    "print(f\"Processing all stations took {time.time() - start_time:.2f} seconds\")\n",
    "\n",
    "# Add the RMSE values to the station_details DataFrame\n",
    "station_details['RMSE_with'] = rmse_values\n",
    "\n",
    "# Save the updated DataFrame to a CSV file\n",
    "station_details.to_csv('rmse_sentinel_highres_with.csv', index=False)\n",
    "\n",
    "# Close the NetCDF file\n",
    "nc_file.close()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "plot",
   "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.18"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}