Spring_verification_codes.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Spring Event Verification"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 1 - Import Necessary Modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cartopy \n",
    "import cartopy.io.shapereader as shpreader \n",
    "import cartopy.io.shapereader as Reader \n",
    "import geopandas as gpd\n",
    "import matplotlib.pyplot as plt \n",
    "import numpy as np \n",
    "import proplot\n",
    "import rioxarray\n",
    "import salem\n",
    "import xarray as xr \n",
    "\n",
    "from cartopy.feature import ShapelyFeature \n",
    "from mpl_toolkits.axes_grid1.inset_locator import inset_axes \n",
    "from netCDF4 import Dataset \n",
    "from shapely.geometry import mapping\n",
    "from wrf import xy_to_ll \n",
    "\n",
    "# warnings\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2 - Station Information"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "station_detail_prec = {   '17130': [39.9727, 32.8637, 'Ankara Bölge'],\n",
    "                     '17134': [39.8032, 32.8434, 'Ufuk Danışment'],\n",
    "                     #'17137': [39.7985, 32.9716, 'Elmadağ Radar'],\n",
    "                     '17651' : [39.1525, 32.1283, 'Polatlı Tigem'],\n",
    "                     '17664' : [40.4729, 32.6441, 'Kızılcahamam'],\n",
    "                     '17679' : [40.1733, 31.332, 'Nallıhan'],\n",
    "                     '17680' : [40.1608, 31.9172, 'Beypazarı'],\n",
    "                     '17715' : [39.92, 33.2125, 'Elmadağ Barutsan Fabrikası'],\n",
    "                     '17728' : [39.5834, 32.1624, 'Polatlı'],\n",
    "                     '17729' : [39.5546, 33.1089, 'Bala'],\n",
    "                     '17731' : [38.9539, 33.4218, 'Şereflikoçhisar'],\n",
    "                     '17733' : [39.613, 32.672, 'Haymana Tarım'],\n",
    "                     '17759' : [39.7414, 32.38, 'Temelli'],\n",
    "                     '17987' : [40.0386, 33.2908, 'Yeşildere'],\n",
    "                     '18045' : [40.03, 32.2345, 'Ayaş'],\n",
    "                     '18046' : [39.9075, 32.8494, 'TBMM Bahçesi'],\n",
    "                     '18072' : [40.0878, 32.8111, 'Bağlum'],\n",
    "                     '18073' : [40.2725, 32.6897, 'Kurtboğazı Barajı'],\n",
    "                     '18074' : [39.9194, 32.9944, 'Bayındır Barajı'],\n",
    "                     '18075' : [40.4839, 32.4608, 'Çamlıdere Orman İşletme Müdürlüğü'],\n",
    "                     '18076' : [40.0947, 33.4133, 'Kalecik DSİ'],\n",
    "                     '18077' : [40.2111, 32.2472, 'Güdül Gençlik Spor'],\n",
    "                     '18078' : [39.5578, 31.9019, 'Acıkır'],\n",
    "                     '18240' : [40.1408, 33.1081, 'Akyurt'],\n",
    "                     '18241' : [40.5983, 32.5033, 'Çamkoru'],\n",
    "                     '18242' : [40.2867, 33.0108, 'Çubuk'],\n",
    "                     '18243' : [40.0317, 32.8933, 'Pursaklar'],\n",
    "                     '18244' : [39.9986, 32.5814, 'Sincan'],\n",
    "                     '18250' : [39.9725, 32.8639, '9. Ankara Bölge'],\n",
    "                     '18256' : [39.8572, 32.8219, 'Çaldağ'],\n",
    "                     '18257' : [39.4356, 32.5153, 'Haymana'],\n",
    "                     '18903' : [40.0761, 32.6075, 'Sarayköy Toprak']}\n",
    "\n",
    "station_detail_temp = {'17127' : [40.0788, 32.5657, 'Akıncı Airport'],\n",
    "                     '17128' : [40.124, 32.9992, 'Esenboğa Airport'],\n",
    "                     '17129' : [39.9558, 32.6854, 'Etimesgut Airport'],\n",
    "                     '17130': [39.9727, 32.8637, 'Ankara Bölge'],\n",
    "                     '17131' : [39.9343, 32.7387, 'Güvercinlik Airport'],\n",
    "                     '17134': [39.8032, 32.8434, 'Ufuk Danışment'],\n",
    "                     '17137': [39.7985, 32.9716, 'Elmadağ Radar'],\n",
    "                     '17651' : [39.1525, 32.1283, 'Polatlı Tigem'],\n",
    "                     '17664' : [40.4729, 32.6441, 'Kızılcahamam'],\n",
    "                     '17679' : [40.1733, 31.332, 'Nallıhan'],\n",
    "                     '17680' : [40.1608, 31.9172, 'Beypazarı'],\n",
    "                     '17715' : [39.92, 33.2125, 'Elmadağ Barutsan Fabrikası'],\n",
    "                     '17728' : [39.5834, 32.1624, 'Polatlı'],\n",
    "                     '17729' : [39.5546, 33.1089, 'Bala'],\n",
    "                     '17731' : [38.9539, 33.4218, 'Şereflikoçhisar'],\n",
    "                     '17733' : [39.613, 32.672, 'Haymana Tarım']}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3 - Open Datasets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#aws obs precipitatio\n",
    "path_prec = r'\\datasets\\ankara_P_tümveri_0918.xlsx'\n",
    "dsp = pd.read_excel(path_prec, sheet_name='ankara_1318_p')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#aws obs temperature\n",
    "path_temp = r'\\datasets\\ankara_0918T.xlsx'\n",
    "dst = pd.read_excel(path_temp, sheet_name='ankara_0918t', na_values=-999)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#SPRING\n",
    "#thompson\n",
    "wrf1_data = salem.open_wrf_dataset(r'\\spring_datasets\\ankara_thompson_urban\\wrfout_d03_2014-05-04_00_00_00') \n",
    "wrf1_data_addit = Dataset(r'\\spring_datasets\\ankara_thompson_urban\\wrfout_d03_2014-05-04_00_00_00')\n",
    "#lin\n",
    "wrf2_data = salem.open_wrf_dataset(r'\\spring_datasets\\ankara_lin_urban\\wrfout_d03_2014-05-04_00_00_00')\n",
    "wrf2_data_addit = Dataset(r'\\spring_datasets\\ankara_lin_urban\\wrfout_d03_2014-05-04_00_00_00')\n",
    "#wsm6\n",
    "wrf3_data = salem.open_wrf_dataset(r'\\spring_datasets\\ankara_wsm6_urban\\wrfout_d03_2014-05-04_00_00_00')\n",
    "wrf3_data_addit = Dataset(r'\\spring_datasets\\ankara_wsm6_urban\\wrfout_d03_2014-05-04_00_00_00')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4 - Extract desired hourly AWS observation Dataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def extract_hourly_point_obs_dataframe(data, start, end):\n",
    "    \"\"\"\n",
    "    data: Pandas Data\n",
    "    start: desired start date (e.g., '2009-01-01 00')\n",
    "    end: desired end date\n",
    "    \"\"\"\n",
    "    data_copy = data.copy(deep=True)\n",
    "    data_start = str(data['year'].iloc[0]) + str(data['month'].iloc[0]) + str(data['day'].iloc[0]) + str(data['hour'].iloc[0])\n",
    "    data_end = str(data['year'].iloc[-1]) + str(data['month'].iloc[-1]) + str(data['day'].iloc[-1]) + str(data['hour'].iloc[-1])\n",
    "    data_start = pd.to_datetime(int(data_start), format='%Y%m%d%H')\n",
    "    data_end = pd.to_datetime(int(data_end), format='%Y%m%d%H')\n",
    "\n",
    "    data_copy.index = pd.date_range(data_start, data_end, freq='1H') #set datetime index\n",
    "    del data_copy['year'], data_copy['month'], data_copy['day'], data_copy['hour']    \n",
    "    data_copy = data_copy.loc[start : end]\n",
    "    return data_copy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#extract point AWS precipitation dataset for the date of interest\n",
    "obs_aws_prec = extract_hourly_point_obs_dataframe(dsp, '2014-05-05 00', '2014-05-05 23')\n",
    "#optional\n",
    "del obs_aws_prec[17725] #Has NaN values ; useless\n",
    "del obs_aws_prec[17137] #Elmadağ Radar outlier in WRF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "#extract point AWS temperature dataset for the date of interest\n",
    "obs_aws_temp = extract_hourly_point_obs_dataframe(dst, '2014-05-05 00', '2014-05-05 23')\n",
    "#if any station contains NaN values you may think to not use that station and delete it\n",
    "#in my case there are no NaN values for the date of my interest"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5 - Extract desired hourly point-wise WRF data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def extract_hourly_point_wrf_dataframe(data, data_addit, station_details, var, start=False, end=False):\n",
    "    \"\"\"\n",
    "    data: WRF dataset opened with Salem \n",
    "    data_addit: WRF dataset opened with NetCDF4\n",
    "    station_details: Dictionary for the details of \n",
    "                     the observation number and lat-lons (e.g., {'17130': [lat,lon]})\n",
    "    var: variable name\n",
    "    \"\"\"\n",
    "    #get the essential information\n",
    "    data_updated = data[var]\n",
    "    data_start = str(data[var].time[0].values)\n",
    "    data_end = str(data[var].time[-1].values)\n",
    "    \n",
    "    #if start and end is given by user\n",
    "    if start and end: \n",
    "        data_start = start\n",
    "        data_end = end\n",
    "        data_updated = data[var].sel(time=slice(start, end)) #update the dataset\n",
    "\n",
    "    data_station_locs = [] #list to create pandas dataFrame\n",
    "    for i in station_details:\n",
    "        x,y = ll_to_xy(data_addit, station_details[i][0], station_details[i][1]) #get x-y values given the station's coordinates\n",
    "        data_station_locs.append(data_updated.isel(south_north=y, west_east=x).values)\n",
    " \n",
    "    data_station_locs = pd.DataFrame(np.array(data_station_locs).transpose(), columns=[int(i) for i in station_details])\n",
    "    data_station_locs.index = pd.date_range(data_start, data_end, freq='1H') #set datetime index\n",
    "    return data_station_locs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Extract WRF point precipitation pandas data for the day of interest\n",
    "wrf1_pandas_prec = extract_hourly_point_wrf_dataframe(wrf1_data, wrf1_data_addit, station_detail_prec, 'PRCP',\n",
    "                                                      start='2014-05-05T00', end='2014-05-05T23')\n",
    "wrf2_pandas_prec = extract_hourly_point_wrf_dataframe(wrf2_data, wrf2_data_addit, station_detail_prec, 'PRCP',\n",
    "                                                      start='2014-05-05T00', end='2014-05-05T23')\n",
    "wrf3_pandas_prec = extract_hourly_point_wrf_dataframe(wrf3_data, wrf3_data_addit, station_detail_prec, 'PRCP',\n",
    "                                                      start='2014-05-05T00', end='2014-05-05T23')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Extract WRF point temperature pandas data for the day of interest\n",
    "wrf1_pandas_temp = extract_hourly_point_wrf_dataframe(wrf1_data, wrf1_data_addit, station_detail_temp, 'T2',\n",
    "                                                      start='2014-05-05T00', end='2014-05-05T23')\n",
    "wrf2_pandas_temp = extract_hourly_point_wrf_dataframe(wrf2_data, wrf2_data_addit, station_detail_temp, 'T2',\n",
    "                                                      start='2014-05-05T00', end='2014-05-05T23')\n",
    "wrf3_pandas_temp = extract_hourly_point_wrf_dataframe(wrf3_data, wrf3_data_addit, station_detail_temp, 'T2',\n",
    "                                                      start='2014-05-05T00', end='2014-05-05T23')\n",
    "\n",
    "wrf1_pandas_temp = wrf1_pandas_temp - 273.15\n",
    "wrf2_pandas_temp = wrf2_pandas_temp - 273.15\n",
    "wrf3_pandas_temp = wrf3_pandas_temp - 273.15"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6 - Calculate Verification Scores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "class metricer(object):\n",
    "    def __init__(self, actual, predicted):\n",
    "        self.actual = actual\n",
    "        self.predicted = predicted\n",
    "    def __repr__(self):\n",
    "        return \"Shape of the actual observation array{} ; Shape of the predicted array{} \".format(np.shape(self.actual), np.shape(self.predicted))\n",
    "    def calc_MAE(self):\n",
    "        \"\"\"\n",
    "        Calculate Mean Absolute Error\n",
    "        \"\"\"\n",
    "        return np.mean(np.abs(self.predicted - self.actual))\n",
    "    \n",
    "    def calc_RMSE(self):\n",
    "        \"\"\"\n",
    "        Calculate Root Mean Squared Error\n",
    "        \"\"\"\n",
    "        return np.sqrt(np.mean(np.square(self.predicted - self.actual)))\n",
    "    \n",
    "    def calc_Pearson_Corr(self):\n",
    "        \"\"\"\n",
    "        Calculate Pearson Correlation Coef\n",
    "        \"\"\"\n",
    "        return np.corrcoef(self.predicted, self.actual)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6.1 - Precipitation Metrics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "#create metric object for each WRF run\n",
    "metric_obj_thompson_prec = metricer(obs_aws_prec, wrf1_pandas_prec)\n",
    "metric_obj_lin_prec = metricer(obs_aws_prec, wrf2_pandas_prec)\n",
    "metric_obj_wsm6_prec = metricer(obs_aws_prec, wrf3_pandas_prec)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.1829473584202836\n",
      "1.5625518077647411\n",
      "1.4222390507547493\n"
     ]
    }
   ],
   "source": [
    "print(metric_obj_thompson_prec.calc_MAE().mean())\n",
    "print(metric_obj_lin_prec.calc_MAE().mean())\n",
    "print(metric_obj_wsm6_prec.calc_MAE().mean())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.8240602300573407\n",
      "2.5510899526808735\n",
      "2.3739742002582336\n"
     ]
    }
   ],
   "source": [
    "print(metric_obj_thompson_prec.calc_RMSE().mean())\n",
    "print(metric_obj_lin_prec.calc_RMSE().mean())\n",
    "print(metric_obj_wsm6_prec.calc_RMSE().mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6.2 - Temperature Metrics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "#create metric object for each WRF run\n",
    "metric_obj_thompson_temp = metricer(obs_aws_temp, wrf1_pandas_temp)\n",
    "metric_obj_lin_temp = metricer(obs_aws_temp, wrf2_pandas_temp)\n",
    "metric_obj_wsm6_temp = metricer(obs_aws_temp, wrf3_pandas_temp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.1157096227010093\n",
      "1.3034344196319578\n",
      "1.5692001660664878\n"
     ]
    }
   ],
   "source": [
    "#Mean Absolute Error MAE\n",
    "print(metric_obj_thompson_temp.calc_MAE().mean())\n",
    "print(metric_obj_lin_temp.calc_MAE().mean())\n",
    "print(metric_obj_wsm6_temp.calc_MAE().mean())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2.368632845647286\n",
      "1.6073513593968216\n",
      "1.8541704568874122\n"
     ]
    }
   ],
   "source": [
    "#Root Mean Square Error RMSE\n",
    "print(metric_obj_thompson_temp.calc_RMSE().mean())\n",
    "print(metric_obj_lin_temp.calc_RMSE().mean())\n",
    "print(metric_obj_wsm6_temp.calc_RMSE().mean())"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}