CR.xrayDS_format.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Aggregating data according to AOD, WV, Precip, and Tau\n",
    "# Formatted data used for CalRad paper\n",
    "####################################################################################################\n",
    "\n",
    "# must insert this statement to render the plots within the notebook\n",
    "# this is specific to the ipython notebook\n",
    "%matplotlib inline\n",
    "\n",
    "## Import statements\n",
    "# netcdf/numpy/xray/stats\n",
    "import numpy as np\n",
    "from datetime import datetime, timedelta\n",
    "import pandas as pd\n",
    "import xarray as xray\n",
    "import kray\n",
    "from scipy import interpolate as interp\n",
    "from scipy.stats.stats import pearsonr\n",
    "from scipy import stats\n",
    "\n",
    "# OS interaction\n",
    "import sys, pickle, os\n",
    "\n",
    "# import subplots function for plotting\n",
    "import seaborn as sns\n",
    "import matplotlib\n",
    "from matplotlib.pyplot import subplots\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import cm\n",
    "from mpl_toolkits.basemap import Basemap\n",
    "import kgraph\n",
    "\n",
    "# Solargeometry\n",
    "import solargeo\n",
    "\n",
    "## Directory listing\n",
    "dir_data = '/Users/karllapo/gdrive/SnowHydrology/proj/CloudClimatology/data'\n",
    "dir_print = '/Users/karllapo/gdrive/SnowHydrology/proj/CloudClimatology/Graphics'\n",
    "\n",
    "# List of sub-directory names for each data set\n",
    "dir_NLDAS = '/NLDAS'\n",
    "dir_SYN = '/CERES_SYN'\n",
    "dir_grobs = '/GroundObs'\n",
    "dir_VIC = '/VIC_MTCLIM'\n",
    "dir_WRF = '/WRF'\n",
    "dir_AN = '/aeronet'\n",
    "\n",
    "# Directory for basemap pickle files\n",
    "dir_bmap = '/Users/karllapo/gdrive/SnowHydrology/proj/CloudClimatology/data/basemap'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "##### Load daily data\n",
    "\n",
    "###########\n",
    "## NLDAS ##\n",
    "os.chdir(dir_data+dir_NLDAS)\n",
    "nldas = xray.open_dataset('CA.NLDAS.irrad.daily.nc')\n",
    "nldas = nldas.rename({'DLWRF_110_SFC':'LWdwn','DSWRF_110_SFC':'SWdwn'})\n",
    "nldas.SWdwn.values[nldas.SWdwn.values > 2000] = np.nan\n",
    "\n",
    "#########\n",
    "## SYN ##\n",
    "os.chdir(dir_data+dir_SYN)\n",
    "syn = xray.open_dataset('CA.SYN.irrad.daily.nc')\n",
    "############ No flipping -- SYN reformatted on 08/20/15 for CR.Composite. OTHER SCRIPTS WILL NEED THIS FIX\n",
    "\n",
    "#########\n",
    "## VIC ##\n",
    "os.chdir(dir_data+dir_VIC)\n",
    "mtclim = xray.open_dataset('CA.MTCLIM.irrad.daily.nc')\n",
    "\n",
    "#########\n",
    "## WRF ##\n",
    "os.chdir(dir_data+dir_WRF)\n",
    "wrf = xray.open_dataset('CA.WRF.irrad.daily.nc')\n",
    "wrf = wrf.rename({'lon':'longitude','lat':'latitude'})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "##### Process daily ground observations\n",
    "os.chdir(dir_data+dir_grobs)\n",
    "grobs = xray.open_dataset('CA.grobs_combined.daily.xray.nc')\n",
    "\n",
    "# Reformat names/coordinates\n",
    "grobs = grobs.resample(freq='D', dim='time', how='mean')\n",
    "grobs.SWdwn.values[grobs.SWdwn.values == 0] = np.nan\n",
    "grobs = grobs.rename({'lon':'longitude','lat':'latitude'})\n",
    "grobs.longitude.values = -grobs.longitude.values\n",
    "\n",
    "# Conversion from python 2.7 to 3.5 caused the xarray strings to become byte literals\n",
    "grobs.network.values = list(map(lambda x: x.decode('utf-8'),grobs.network.values))\n",
    "grobs.station.values = list(map(lambda x: x.decode('utf-8'),grobs.station.values))\n",
    "\n",
    "## Remove stations outside study domain\n",
    "# Bounding box - ragged domain\n",
    "LL_rag = [-120,34.5]\n",
    "LR_rag = [-115,34.5]\n",
    "UR_rag = [-118.5,41]\n",
    "UL_rag = [-123.5,41]\n",
    "\n",
    "# Station lat/lon\n",
    "stat_lat = grobs.latitude.values\n",
    "stat_lon = grobs.longitude.values\n",
    "\n",
    "# Ragged domain, CA study area\n",
    "line_west_m = (UL_rag[1]-LL_rag[1])/(UL_rag[0]-LL_rag[0])\n",
    "line_west_b = LL_rag[1]-line_west_m*LL_rag[0]\n",
    "line_east_m = (UR_rag[1]-LR_rag[1])/(UR_rag[0]-LR_rag[0])\n",
    "line_east_b = LR_rag[1]-line_east_m*LR_rag[0]\n",
    "ind = np.nonzero((stat_lon > (stat_lat -line_west_b)/line_west_m) & \\\n",
    "                (stat_lon < (stat_lat-line_east_b)/line_east_m) & \\\n",
    "                (stat_lat > LR_rag[1]) & (stat_lat < UL_rag[1]))\n",
    "# Reindex\n",
    "stat_to_keep = grobs.station[ind]\n",
    "grobs = grobs.reindex(station=stat_to_keep)\n",
    "\n",
    "## Append station group information\n",
    "# Station data excel speadsheet\n",
    "stdat = pd.read_csv('All_StationSummary.v2.csv',sep= ',', \\\n",
    "                index_col=0,na_values=[-9999,'NaN']) # Read the supporting information\n",
    "# station data -> xray structure\n",
    "stdat = xray.Dataset(coords={'station': (['station'], stdat.index), \\\n",
    "                             'Grouping': (['station'],  stdat.Grouping)})    \n",
    "# Merge into ground observation structure\n",
    "grobs = grobs.merge(stdat,join='inner')\n",
    "grobs = grobs.rename({'SWdwn':'grobs'})\n",
    "\n",
    "# Removed bad stations, but netcdf no longer contains these stations\n",
    "grobs = grobs.drop(['smj','Pacific_Grove','Woodland'],dim='station')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#### Nan months with less than 90% of days observing\n",
    "numdays_permonth = (~np.isnan(grobs.grobs)).resample(freq='M', dim='time', how='mean',label='right')\n",
    "numdays_permonth=numdays_permonth.reindex(time=grobs.time,method='bfill')\n",
    "grobs.grobs.values[numdays_permonth.values < .9] = np.nan"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#### Combine\n",
    "## List w/ all irradiance datasets\n",
    "daily_mean = {}\n",
    "daily_mean['syn'] = syn\n",
    "daily_mean['nldas'] = nldas\n",
    "daily_mean['mtclim'] = mtclim\n",
    "daily_mean['wrf'] = wrf\n",
    "daily_mean['grobs'] = grobs\n",
    "\n",
    "## Names\n",
    "# Product names\n",
    "pr_names = ['mtclim','nldas','syn','wrf']\n",
    "num_products = np.size(pr_names)\n",
    "grouping = ['north cv','south cv','west of crest','foothills','east of crest']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "####################################################\n",
    "## Find grid point containing each ground station ##\n",
    "####################################################\n",
    "# Station lat and lon\n",
    "lon_stat = grobs.longitude.values\n",
    "lat_stat = grobs.latitude.values\n",
    "\n",
    "for pr in pr_names:        \n",
    "    if pr == 'grobs':\n",
    "        continue\n",
    "    # lat/lon for product\n",
    "    lon_rad = daily_mean[pr].longitude.values\n",
    "    lat_rad = daily_mean[pr].latitude.values  \n",
    "    # mesh\n",
    "    lonm, latm = np.meshgrid(lon_rad,lat_rad)\n",
    "    \n",
    "    # Empty numpy array\n",
    "    to_merge = np.empty((daily_mean[pr].time.size,grobs.station.size))\n",
    "    \n",
    "    ## Product values in each grid containing station\n",
    "    for stat in grobs.station.values:\n",
    "        \n",
    "        # Station index\n",
    "        stat_ind = np.where(stat == grobs.station.values)\n",
    "        # Distance to product grid lat-lon\n",
    "        d = (latm-lat_stat[stat_ind])**2 + (lonm-lon_stat[stat_ind])**2\n",
    "        # Index of closest product grid\n",
    "        dind = np.where(d==np.amin(d))\n",
    "        # Grad grid values at the station, put into xray dataset\n",
    "        if (daily_mean[pr].SWdwn.values[:,dind[0][0],dind[1][0]] == 0).all() \\\n",
    "                | np.isnan(daily_mean[pr].SWdwn.values[:,dind[0][0],dind[1][0]]).all():\n",
    "            to_merge[:,stat_ind[0]] = daily_mean[pr].SWdwn.values[:,dind[0][0]+1,dind[1][0],np.newaxis]\n",
    "        else:\n",
    "            to_merge[:,stat_ind[0]] = daily_mean[pr].SWdwn.values[:,dind[0][0],dind[1][0],np.newaxis]\n",
    "    \n",
    "    ## Merge products w/ grobs xray structure\n",
    "    to_merge_ds = xray.Dataset({pr:(('time','station'),to_merge), \\\n",
    "                                    'time':daily_mean[pr].time.values,\\\n",
    "                                    'station':grobs.station.values})\n",
    "    grobs = grobs.merge(to_merge_ds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#### Mean bias for each group\n",
    "grobs_diff = xray.Dataset()\n",
    "grouping_var = 'Grouping'\n",
    "for pr in pr_names:\n",
    "    grdiff = grobs[pr]-grobs['grobs']\n",
    "    ## Mean bias for each group\n",
    "    grdiff = kray.group_mean(grdiff.swap_dims({'station':'Grouping'}),grouping_var,grouping)\n",
    "    grobs_diff[pr] = kray.combinevars(grdiff,grouping,new_dim_name='Grouping',combinevarname=pr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "##   Transmissivity (products and ground obs.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "//anaconda/lib/python3.5/site-packages/ipykernel/__main__.py:5: FutureWarning: how in .resample() is deprecated\n",
      "the new syntax is .resample(...).mean()\n"
     ]
    }
   ],
   "source": [
    "#### Elevation Angle\n",
    "d = pd.to_datetime(grobs.time.values)\n",
    "el = xray.Dataset()\n",
    "for stat in grobs.station.values:\n",
    "    el[stat] = (('time'),solargeo.avg_el(d,\\\n",
    "                     grobs.loc[{'station':stat}].latitude.values,\\\n",
    "                     grobs.loc[{'station':stat}].longitude.values,\\\n",
    "                     ref='END'))\n",
    "el = kray.combinevars(el,el.data_vars,new_dim_name='station',combinevarname='el')\n",
    "el.coords['time'] = d\n",
    "\n",
    "# Add to xray Dataset\n",
    "grobs['el'] = (('station','time'),el)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Transmissivity\n",
    "tau = grobs.grobs/(np.sin(grobs.el*np.pi/180)*1365)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#### Monthly anomaly\n",
    "ds_in = tau.swap_dims({'station':'Grouping'})\n",
    "grouping_var = 'Grouping'\n",
    "var = 'tau'\n",
    "period = 'month'\n",
    "\n",
    "anom_mon = kray.group_anom(ds_in,grouping_var,grouping,var,period)\n",
    "anom_mon = anom_mon.resample(freq='M',dim='time',how='mean')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#### Daily anomaly\n",
    "ds_in = tau.swap_dims({'station':'Grouping'})\n",
    "grouping_var = 'Grouping'\n",
    "var = 'tau'\n",
    "period = 'day'\n",
    "\n",
    "anom_day = kray.group_anom(ds_in,grouping_var,grouping,var,period)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "##     Auxilary data @ Ground Observations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#### Precip obs @ CIMIS stations\n",
    "# Load\n",
    "os.chdir(dir_data+dir_grobs)\n",
    "precip = xray.open_dataset('CA.CIMIS.Precip.daily.xray.nc')\n",
    "\n",
    "# Conversion from python 2.7 to 3.5 caused the xarray strings to become byte literals\n",
    "precip.network.values = list(map(lambda x: x.decode('utf-8'),precip.network.values))\n",
    "precip.station.values = list(map(lambda x: x.decode('utf-8'),precip.station.values))\n",
    "\n",
    "# Format\n",
    "precip = precip.Precip\n",
    "precip = precip.resample(freq='D', dim='time', how='mean')\n",
    "\n",
    "# Add to xray Dataset\n",
    "precip = precip.reindex_like(grobs)\n",
    "grobs['precip'] = (('station','time'),precip)\n",
    "\n",
    "#### Aggregate to grouping variables\n",
    "ds_in = grobs.precip.swap_dims({'station':'Grouping'})\n",
    "grouping_var = 'Grouping'\n",
    "\n",
    "precip = kray.group_mean(ds_in,grouping_var,grouping)\n",
    "precip = kray.combinevars(precip,grouping,new_dim_name='Grouping',combinevarname='precip')\n",
    "\n",
    "## Add to anom xray data structure\n",
    "anom_day['precip'] = (('Grouping','time'),precip)\n",
    "grobs_diff['precip'] = (('Grouping','time'),precip)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "##### total aerosol optical depth\n",
    "def parse(str1, str2):\n",
    "    date_str = str1+'_'+str2\n",
    "    dt = datetime.strptime(date_str,\"%d:%m:%Y_%H:%M:%S\")\n",
    "    return dt\n",
    "\n",
    "iter_count = 0\n",
    "os.chdir(dir_data+dir_AN)\n",
    "fid = '020101_121231_Fresno.ONEILL_20'\n",
    "with open(fid, 'r') as datafile:\n",
    "    # Skip the header of arbitrary size and read the column names\n",
    "    if iter_count == 0:\n",
    "        line = datafile.readline()\n",
    "        iter_count = iter_count+1\n",
    "    while not line.startswith('Date(dd:mm:yyyy)'):\n",
    "        line = datafile.readline()\n",
    "    line = line.replace('\\n','')\n",
    "    col_names = line.split(',')\n",
    "    col_names[-1] = 'wavelength-1'\n",
    "    col_names.append('wavelength-2')\n",
    "    col_names.append('wavelength-3')\n",
    "    col_names.append('wavelength-4')\n",
    "    col_names.append('wavelength-5')\n",
    "    data = pd.read_csv(datafile, names=col_names, sep= ',',\\\n",
    "                       parse_dates={'Datetime' : [0,1]},date_parser=parse,\\\n",
    "                       index_col='Datetime')\n",
    "aod = data['Total_AOD_500nm[tau_a]']\n",
    "\n",
    "## Add to xray Dataset\n",
    "aod = aod.reindex(grobs.time.values,fill_value=np.nan)\n",
    "grobs['aod'] = (('time'),aod)\n",
    "anom_day['aod'] = (('time'),aod)\n",
    "grobs_diff['aod'] = (('time'),aod)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "##### water vapor\n",
    "iter_count = 0\n",
    "fid = '020101_121231_Fresno.lev20'\n",
    "with open(fid, 'r') as datafile:\n",
    "    # Skip the header of arbitrary size and read the column names\n",
    "    if iter_count == 0:\n",
    "        line = datafile.readline()\n",
    "        iter_count = iter_count+1\n",
    "    while not line.startswith('Date(dd-mm-yy)'):\n",
    "        iter_count = iter_count+1\n",
    "        line = datafile.readline()\n",
    "        if iter_count > 10:\n",
    "            break\n",
    "    line = line.replace('\\n','')\n",
    "    col_names = line.split(',')\n",
    "    data = pd.read_csv(datafile, names=col_names, sep= ',',\\\n",
    "                       parse_dates={'Datetime' : [0,1]},date_parser=parse,\\\n",
    "                       index_col='Datetime')\n",
    "wv = data['Water(cm)']\n",
    "\n",
    "## Add to xray Dataset\n",
    "wv = wv.reindex(grobs.time.values,fill_value=np.nan)\n",
    "grobs['wv'] = (('time'),wv)\n",
    "anom_day['wv'] = (('time'),wv)\n",
    "grobs_diff['wv'] = (('time'),wv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#### Save formatted data for later use in other scripts\n",
    "os.chdir(dir_data)\n",
    "grobs.to_netcdf('grobs.daily.nc',\n",
    "                encoding={'station': {'dtype': 'str'},\n",
    "                          'network': {'dtype': 'str'}})\n",
    "grobs_diff.to_netcdf('grobs_bias.daily.nc')\n",
    "anom_day.to_netcdf('tau_anom.daily.nc')\n",
    "anom_mon.to_netcdf('tau_anom.monthly.nc')\n",
    "tau.to_dataset(name='tau').to_netcdf('tau.daily.nc',\n",
    "                                     encoding={'station': {'dtype': 'str'},\n",
    "                                               'network': {'dtype': 'str'}})"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [Root]",
   "language": "python",
   "name": "Python [Root]"
  },
  "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}