PVSC44 TL sensitivity.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "import h5py\n",
    "from matplotlib import pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import pvlib\n",
    "import seaborn as sns\n",
    "import statsmodels.api as sm\n",
    "\n",
    "from pvsc44_clearsky_aod import ecmwf_macc_tools\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "sns.set_context('notebook', rc={'figure.figsize': (8, 6)})\n",
    "sns.set(font_scale=1.5)\n",
    "\n",
    "# get the \"metadata\" that contains the station id codes for the SURFRAD data that was analyzed\n",
    "METADATA = pd.read_csv('metadata.csv', index_col=0)\n",
    "\n",
    "\n",
    "# load calculations for each station\n",
    "atm_params_3min_clear = {}\n",
    "for station_id in METADATA.index:\n",
    "    with h5py.File('%s_3min_clear_atm_params.h5' % station_id, 'r') as f:\n",
    "        np_atm_params_3min_clear = pd.DataFrame(np.array(f['data']))\n",
    "    np_atm_params_3min_clear['index'] = pd.DatetimeIndex(np_atm_params_3min_clear['index'])\n",
    "    np_atm_params_3min_clear.set_index('index', inplace=True)\n",
    "    np_atm_params_3min_clear.index.rename('timestamps', inplace=True)    \n",
    "    atm_params_3min_clear[station_id] = np_atm_params_3min_clear\n",
    "    \n",
    "    \n",
    "    \n",
    "# filter out low light\n",
    "\n",
    "# CONSTANTS\n",
    "MODELS = {'solis': 'SOLIS', 'lt': 'Linke', 'macc': 'ECMWF-MACC', 'bird': 'Bird'}\n",
    "CS = ['dni', 'dhi', 'ghi']\n",
    "LOW_LIGHT = 200  # threshold for low light in W/m^2\n",
    "\n",
    "is_bright = {}\n",
    "for station_id, station_atm_params_3min_clear in atm_params_3min_clear.iteritems():\n",
    "    is_bright[station_id] = station_atm_params_3min_clear['ghi'] > LOW_LIGHT\n",
    "    \n",
    "    \n",
    "    \n",
    "TIMES = pd.DatetimeIndex(start='2003-01-01T00:00:00', freq='3H', end='2012-12-31T23:59:59').tz_localize('UTC')\n",
    "LT = pvlib.clearsky.lookup_linke_turbidity(\n",
    "    time=TIMES,\n",
    "    latitude=METADATA['latitude'][station_id],\n",
    "    longitude=METADATA['longitude'][station_id],)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "            utc offset           tz    station name  latitude  longitude  \\\n",
      "station id                                                                 \n",
      "bon               -6.0   US/Central       Bondville     40.05     -88.37   \n",
      "tbl               -7.0  US/Mountain  Table Mountain     40.13    -105.24   \n",
      "dra               -8.0   US/Pacific     Desert Rock     36.62    -116.02   \n",
      "fpk               -7.0  US/Mountain       Fort Peck     48.31    -105.10   \n",
      "gwn               -6.0   US/Central   Goodwin Creek     34.25     -89.87   \n",
      "psu               -5.0   US/Eastern      Penn State     40.72     -77.93   \n",
      "sxf               -6.0   US/Central     Sioux Falls     43.73     -96.62   \n",
      "\n",
      "            elevation  \n",
      "station id             \n",
      "bon               213  \n",
      "tbl              1689  \n",
      "dra              1007  \n",
      "fpk               634  \n",
      "gwn                98  \n",
      "psu               376  \n",
      "sxf               473  \n"
     ]
    }
   ],
   "source": [
    "print METADATA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def change_lt(df_in, value):   \n",
    "    \n",
    "    df = df_in.copy(deep=True)\n",
    "    df = df[df['ghi']>200]\n",
    "    \n",
    "    df['lt_ref'] = df['lt']\n",
    "    df['lt'] = value       \n",
    "    # solarposition, airmass, pressure-corrected & water vapor\n",
    "    SOLPOS = pvlib.solarposition.get_solarposition(\n",
    "        time=df.index,\n",
    "        latitude=METADATA['latitude'][station_id],\n",
    "        longitude=METADATA['longitude'][station_id],\n",
    "        altitude=METADATA['elevation'][station_id],\n",
    "        pressure=df['press'],\n",
    "        temperature=df['ta']\n",
    "    )\n",
    "    AIRMASS = pvlib.atmosphere.relativeairmass(SOLPOS.apparent_zenith)  # relative airmass\n",
    "    AM_PRESS = pvlib.atmosphere.absoluteairmass(AIRMASS, pressure=df['press'])  # pressure corrected airmass\n",
    "    ETR = pvlib.irradiance.extraradiation(df.index)  # extra-terrestrial total solar radiation\n",
    "\n",
    "    # append these to the SOLPOS dataframe to keep them together more easily\n",
    "    df['apparent_zenith'] = SOLPOS.apparent_zenith\n",
    "    df['am']= AIRMASS\n",
    "    df['amp']= AM_PRESS\n",
    "\n",
    "    df['etr']= ETR\n",
    "    \n",
    "    \n",
    "    INEICHEN_LT = pvlib.clearsky.ineichen(\n",
    "    df['apparent_zenith'], df['amp'], df['lt'], altitude=METADATA['elevation'][station_id],\n",
    "    dni_extra=df['etr']\n",
    "    )\n",
    "\n",
    "    df['lt_dhi']= INEICHEN_LT['dhi']\n",
    "    df['lt_dni']= INEICHEN_LT['dni']\n",
    "    df['lt_ghi']= INEICHEN_LT['ghi']\n",
    "    \n",
    "    \n",
    "    return df\n",
    "\n",
    "\n",
    "def get_monthly_mbe(df):\n",
    "    \n",
    "    model = 'lt'\n",
    "    \n",
    "    output = pd.DataFrame()\n",
    "    \n",
    "    for cs in CS:\n",
    "        monthly_avg = df[cs].resample('M').mean()\n",
    "        output['mbe_'+ cs] = (df['%s_%s' % (model, cs)] - df[cs]).resample('M').mean() / monthly_avg    \n",
    "        \n",
    "    output['month'] = output.index.month\n",
    "    \n",
    "    return output\n",
    "\n",
    "\n",
    "def find_best_lt(df_in, month, lt_start):\n",
    "    \n",
    "    df = df_in.copy(deep=True)\n",
    "    df = df[df.index.month==month]\n",
    "    \n",
    "    \n",
    "    median_error =10\n",
    "    lt = lt_start\n",
    "    amount = 1.0\n",
    "    \n",
    "    \n",
    "    df = change_lt(df, lt_start)\n",
    "    mbe = get_monthly_mbe(df)\n",
    "    median_error = mbe['mbe_ghi'].median()\n",
    "    \n",
    "    direction = median_error/abs(median_error)\n",
    "    last_direction = median_error/abs(median_error)\n",
    "    \n",
    "    \n",
    "    while abs(median_error) > 0.01:\n",
    "        \n",
    "        direction = median_error/abs(median_error)\n",
    "        \n",
    "        if last_direction != direction:\n",
    "            amount *= 0.5\n",
    "            \n",
    "        last_direction = direction\n",
    "                \n",
    "        \n",
    "        df = change_lt(df, lt)\n",
    "        mbe = get_monthly_mbe(df)\n",
    "        median_error = mbe['mbe_ghi'].median()\n",
    "        \n",
    "        \n",
    "        \n",
    "        \n",
    "        if median_error > 0:\n",
    "            lt += amount\n",
    "        else:\n",
    "            lt -= amount\n",
    "            \n",
    "        #rint direction, median_error, lt, amount\n",
    "            \n",
    "    \n",
    "    #print i\n",
    "    \n",
    "    return lt\n",
    "                \n",
    "        \n",
    "        \n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#make a table of \n",
    "    \n",
    "    \n",
    "def get_lt(station_id):\n",
    "    \n",
    "    \n",
    "    TIMES = pd.DatetimeIndex(start='2003-01-01T00:00:00', freq='3H', end='2012-12-31T23:59:59').tz_localize('UTC')\n",
    "    LT = pvlib.clearsky.lookup_linke_turbidity(\n",
    "        time=TIMES,\n",
    "        latitude=METADATA['latitude'][station_id],\n",
    "        longitude=METADATA['longitude'][station_id],)\n",
    "    lt_df = pd.DataFrame(index=LT.index, data=LT, columns=['lt'])\n",
    "    lt_df['month'] = lt_df.index.month\n",
    "    lt_m = lt_df.groupby(\"month\").median()\n",
    "    return lt_m\n",
    "    \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "bon\n",
      "bon 1 lt_ref 2.35508474576 lt_opt 1.91758474576\n",
      "bon 2 lt_ref 2.655 lt_opt 2.155\n",
      "bon 3 lt_ref 2.95573770492 lt_opt 2.58073770492\n",
      "bon 4 lt_ref 3.31803278689 lt_opt 2.56803278689\n",
      "bon 5 lt_ref 4.35245901639 lt_opt 4.35245901639\n",
      "bon 6 lt_ref 4.29672131148 lt_opt 3.29672131148\n",
      "bon 7 lt_ref 4.12016129032 lt_opt 3.87016129032\n",
      "bon 8 lt_ref 4.13770491803 lt_opt 3.38770491803\n",
      "bon 9 lt_ref 4.08360655738 lt_opt 4.08360655738\n",
      "bon 10 lt_ref 3.09262295082 lt_opt 2.84262295082\n",
      "bon 11 lt_ref 2.64508196721 lt_opt 2.64508196721\n",
      "bon 12 lt_ref 2.35 lt_opt 2.225\n",
      "tbl\n",
      "tbl 1 lt_ref 2.75338983051 lt_opt 0.253389830508\n",
      "tbl 2 lt_ref 2.95166666667 lt_opt 0.951666666667\n",
      "tbl 3 lt_ref 3.05983606557 lt_opt 0.309836065574\n",
      "tbl 4 lt_ref 3.65327868852 lt_opt 2.40327868852\n",
      "tbl 5 lt_ref 3.85327868852 lt_opt 3.85327868852\n",
      "tbl 6 lt_ref 4.05491803279 lt_opt 4.80491803279\n",
      "tbl 7 lt_ref 4.26612903226 lt_opt 4.26612903226\n",
      "tbl 8 lt_ref 3.94426229508 lt_opt 3.94426229508\n",
      "tbl 9 lt_ref 3.59508196721 lt_opt 2.22008196721\n",
      "tbl 10 lt_ref 3.29262295082 lt_opt 1.04262295082\n",
      "tbl 11 lt_ref 2.84836065574 lt_opt 0.348360655738\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\gkimball\\Envs\\env_ecmwf\\lib\\site-packages\\pvlib\\clearsky.py:138: RuntimeWarning: invalid value encountered in minimum\n",
      "  dni = np.minimum(bnci, bnci_2)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "tbl 12 lt_ref 2.75 lt_opt -0.25\n",
      "dra\n",
      "dra 1 lt_ref 2.40338983051 lt_opt -1.15911016949\n"
     ]
    }
   ],
   "source": [
    "\n",
    "\n",
    "summary=pd.DataFrame()\n",
    "\n",
    "\n",
    "for station in METADATA.index:\n",
    "    \n",
    "    print station\n",
    "\n",
    "    for month in range(1,13):\n",
    "        \n",
    "        lt_m = get_lt(station)\n",
    "        lt_ref = lt_m.loc[month,'lt']\n",
    "        \n",
    "        lt_opt = find_best_lt(atm_params_3min_clear[station], month, lt_ref)  \n",
    "        \n",
    "        df = pd.DataFrame(index=[month])\n",
    "        df['station'] = station\n",
    "        df['lt_opt'] = lt_opt\n",
    "        df['lt_ref'] = lt_ref\n",
    "        df['month'] = month\n",
    "\n",
    "        \n",
    "        print station, month, 'lt_ref', lt_ref, 'lt_opt', lt_opt\n",
    "        summary = summary.append(df, ignore_index=True)    \n",
    "        summary.to_csv('summary.csv')\n",
    "\n",
    "    \n",
    "    \n",
    "\n",
    "    \n",
    "    \n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}