RMSEv2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Error retrieving data for station 8311030: CO-OPS API returned an error. Status Code: 400. Reason: Bad Request\n",
      "\n",
      "Error retrieving data for station 8311062: CO-OPS API returned an error. Status Code: 400. Reason: Bad Request\n",
      "\n"
     ]
    }
   ],
   "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",
    "\n",
    "# Function to calculate RMSE\n",
    "def calculate_rmse(true_values, predicted_values):\n",
    "    if len(true_values) > 48 and len(predicted_values) > 48:\n",
    "        true_values = true_values[48:]\n",
    "        predicted_values = predicted_values[48:]\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",
    "station_details = pd.read_csv('station_details.csv')\n",
    "\n",
    "# Load the NetCDF file\n",
    "nc_file = netCDF4.Dataset('/work2/07174/soelem/hopper/fort.63.nc', 'r')\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",
    "kdtree = cKDTree(points)\n",
    "\n",
    "rmse_values = []\n",
    "\n",
    "# Chunk size for reading 'zeta' variable\n",
    "chunk_size = 1000\n",
    "\n",
    "for index, row in station_details.iterrows():\n",
    "    try:\n",
    "        station_id = row['Station ID']\n",
    "        _, closest_point_idx = kdtree.query([row['Latitude'], row['Longitude']])\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",
    "        # Reading 'zeta' data in chunks\n",
    "        nc_zeta_adjusted = []\n",
    "        for i in range(0, len(nc_file.variables['zeta'][:, closest_point_idx]), chunk_size):\n",
    "            chunk = nc_file.variables['zeta'][i:i + chunk_size, closest_point_idx]\n",
    "            chunk = np.delete(chunk, np.s_[0:4]) if i == 0 else chunk  # Delete the first 4 values only for the first chunk\n",
    "            nc_zeta_adjusted.append(chunk)\n",
    "\n",
    "        nc_zeta_adjusted = np.concatenate(nc_zeta_adjusted)\n",
    "\n",
    "        nc_zeta_series = pd.Series(nc_zeta_adjusted.flatten(), index=pd.date_range(start=\"2023-07-16 00:00:00\", periods=len(nc_zeta_adjusted), 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",
    "station_details['RMSE_without'] = rmse_values\n",
    "station_details.to_csv('rmse_jason_highres_without.csv', index=False)\n",
    "\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
}