calc_dfm_gridcell_is.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import xray\n",
    "# from xray import ufuncs\n",
    "import os\n",
    "import math \n",
    "import gc \n",
    "import pandas as pd \n",
    "from snowpack_functions import get_elev_for_lat_lon,import_gridcell_elevation\n",
    "import datetime as dt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "IMPORT DATA FROM INTEGRATED SCENARIOS ARCHIVE "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "7535"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "################################# INPUTS #############################\n",
    "#lat = 48.03125\n",
    "lat = 48.021\n",
    "#lon = -121.28125\n",
    "lon = -121.25\n",
    "######################################################################\n",
    "\n",
    "direc = '/raid/gergel/%s' % \"tmin\"\n",
    "tmin_file = \"%s_%s_%s.nc\" % (model,scenario,\"tasmin\")\n",
    "tmin_full = xray.open_dataset(os.path.join(direc,tmin_file)) ## load tmin\n",
    "\n",
    "direc = '/raid/gergel/%s' % \"tmax\"\n",
    "tmax_file = \"%s_%s_%s.nc\" % (model,scenario,\"tasmax\")\n",
    "tmax_full = xray.open_dataset(os.path.join(direc,tmax_file)) ## load tmax\n",
    "\n",
    "direc = '/raid/gergel/%s' % \"rh\"\n",
    "q_file = \"%s_%s_%s.nc\" % (model,scenario,\"huss\")\n",
    "q_full = xray.open_dataset(os.path.join(direc,rhmax_file)) ## load specific humidity\n",
    "\n",
    "direc = '/raid/gergel/%s' % \"precip\"\n",
    "pr_file = \"%s_%s_%s.nc\" % (model,scenario,\"pr\")\n",
    "pr_full = xray.open_dataset(os.path.join(direc,pr_file)) ## load precip\n",
    "\n",
    "## adjust lat/lon dimensions since the index names are different\n",
    "tmin_lons_new = tmin_full['lon'].values[tmin_full['lon'].values > 180] - 360 \n",
    "tmin_full['lon'] = tmin_lons_new\n",
    "tmax_full['lon'] = tmin_lons_new\n",
    "pr_full['lon'] = tmin_lons_new\n",
    "q_full['lon'] = tmin_lons_new \n",
    "\n",
    "## select out gridcell \n",
    "lon_ind = np.argmin(np.abs(pr_full.lon - lon))\n",
    "lat_ind = np.argmin(np.abs(pr_full.lat - lat))\n",
    "\n",
    "## get julian days \n",
    "julians = pd.DatetimeIndex(np.asarray(pr_full.day)).dayofyear\n",
    "\n",
    "tmin = tmin_full.isel_points(lon=[lon_ind],lat=[lat_ind])['air_temp_min'].values[0,:] ## air_temp_min\n",
    "tmax = tmax_full.isel_points(lon=[lon_ind],lat=[lat_ind])['air_temp_max'].values[0,:] ## air_temp_max\n",
    "q = q_full.isel_points(lon=[lon_ind],lat=[lat_ind])['specific_humidity'].values[0,:] ## specific_humidity\n",
    "pptamt = pr_full.isel_points(lon=[lon_ind],lat=[lat_ind]) ## variable is precipitation \n",
    "\n",
    "del tmax_full, q_full,pr_full,tmin_full\n",
    "gc.collect() "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "100-hr and 1000-hr DFM FUNCTION "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def calc_fm100_fm1000(pptdur,maxrh,minrh,maxt,mint,lat,tmois,bv,julians,ymc100): \n",
    "    \"\"\"this subroutine computes the average boundary conditions for the past \n",
    "    24 hour and 100-hr-tl fuel moisture. The boundary conditions are weighted averages \n",
    "    of the EQMCs calculated from the temp and RH values. Philab is used to calculate \n",
    "    daylength which is the basis of the weighting function.\"\"\" \n",
    "\n",
    "    emc1 = 0\n",
    "    emc2 = 0\n",
    "    emc = 0\n",
    "    bndry1 = 0\n",
    "    bndry = 0\n",
    "    bvave = 0\n",
    "    daylit = 0\n",
    "    ambvp = 0\n",
    "    fr100 = 0.3156\n",
    "\n",
    "    '''# daylit = calcDaylight(j_date,lat) ########## need to write this function \n",
    "    #######################\n",
    "    phi = lat * 0.01745 ## lat is in degrees \n",
    "    decl = 0.41008*np.sin(np.deg2rad((julians-82) * 0.01745))\n",
    "    daylit = 24 * (1.0 - ((np.arccos(np.tan(np.deg2rad(phi)) * np.tan(np.deg2rad(decl))) / math.pi)))'''\n",
    "    #######################\n",
    "    \n",
    "    ## John's calcDaylight function\n",
    "    if julians > 365:\n",
    "        julians = 365\n",
    "    phi = lat * 0.01745 ## converts latitude to radians\n",
    "    decl = .41008*np.sin((julians-82)*0.01745)\n",
    "    daylit = 24.0*(1-np.arccos(np.tan(phi)*np.tan(decl))/3.14159)\n",
    "    daylit = daylit.real\n",
    "    \n",
    "    #emc1 = np.ndarray(shape=(len(maxrh)),dtype='float')\n",
    "    if minrh <= 10:\n",
    "        emc1= 0.03229 + (0.281073 * minrh) - (0.000578 * minrh * maxt)\n",
    "    elif minrh > 10 and minrh <= 50: \n",
    "        emc1 = 2.22749 + (0.160107 * minrh) - (0.014784 * maxt)\n",
    "    else: \n",
    "        emc1 = 21.0606 + (0.005565 * (minrh**2)) - (0.00035 * minrh * maxt) - (0.483199 * minrh)\n",
    "\n",
    "    #emc2 = np.ndarray(shape=(len(maxrh)),dtype='float')\n",
    "    if maxrh <= 10:\n",
    "        emc2 = 0.03229 + (0.281073 * maxrh) - (0.000578 * maxrh* mint)\n",
    "    elif maxrh > 10 and maxrh <= 50: \n",
    "        emc2 = 2.22749 + (0.160107 * maxrh) - (0.014784 * mint)\n",
    "    else: \n",
    "        emc2 = 21.0606 + (0.005565 * (maxrh**2)) - (0.00035 * maxrh * mint) - (0.483199 * maxrh)\n",
    "    \n",
    "    ## qc maxrh \n",
    "    # for day in np.arange(len(emc1)):\n",
    "    if np.isnan(maxrh):\n",
    "        emc1 = np.nan\n",
    "        emc2 = np.nan \n",
    "\n",
    "    emc = (daylit * emc1 + (24.0 - daylit) * emc2) / 24.0 \n",
    "\n",
    "    ## qc precip duration \n",
    "    # for day in np.arange(len(pptdur)):\n",
    "    \n",
    "    # pptdur = pptdur*1.25\n",
    "    \n",
    "    if pptdur < 0:\n",
    "        pptdur = 0\n",
    "    elif pptdur > 8:\n",
    "        pptdur = 8\n",
    "        \n",
    "    bndry1 = ((24.0 - pptdur) * emc + (0.5 * pptdur + 41) * pptdur) / 24.0 \n",
    "    fm100 = ((bndry1 - ymc100) * fr100) + ymc100 \n",
    "    ## calculate 1000-hr fuel moisture daily using average of boundary conditions for\n",
    "    ## past seven days. starting value set by climate type. \n",
    "\n",
    "    fr1 = 0.3068\n",
    "\n",
    "    bvave = 0\n",
    "\n",
    "    ## accumulate a 6-day total\n",
    "    for i in np.arange(0,6):\n",
    "        bv[:,i] = bv[:,i+1]\n",
    "        bvave = bvave + bv[:,i]\n",
    "\n",
    "    bndry = ((24 - pptdur) * emc + (2.7 * pptdur + 76) * pptdur) / 24.0 \n",
    "    bv[:,6] = bndry\n",
    "\n",
    "    ## add today's boundary from subfm100, divide by 7 days \n",
    "    bvave = (bvave + bndry) / 7.0 \n",
    "\n",
    "    ## calculate today's 1000 hr fuel moisture \n",
    "    # fm1000 = tmois[:,1] + (bvave - tmois[:,1])*fr1 \n",
    "    fm1000 = tmois[:,0] + (bvave - tmois[:,0])*fr1\n",
    "\n",
    "    ## move each days 1000 hr down one, drop the oldest \n",
    "    for i in np.arange(0,6):\n",
    "        tmois[:,i] = tmois[:,i+1] \n",
    "    tmois[:,6] = fm1000 \n",
    "\n",
    "    return(tmois,fm1000,fm100,bv)\n",
    "\n",
    "def kelvin_to_fahrenheit(T):\n",
    "    ''' converts T in Kelvin to Fahrenheit'''\n",
    "    F = ((9.0/5.0) * (T - 273.15)) + 32.0\n",
    "    return(F)\n",
    "\n",
    "def estimate_p(h):\n",
    "    p = 101325 * (1 - ((2.25577 * 10**-5) * h))**5.25588 ## Pascals \n",
    "    return(p)\n",
    "\n",
    "def estimate__e_s(T): \n",
    "    ''' estimates saturation vapor pressure'''\n",
    "    T0 = 273.15 ## Kelvin, reference temperature\n",
    "    e_s = 611 * np.exp((17.67 * (T - T0)) / (T - 29.65) )\n",
    "    return(e_s)\n",
    "\n",
    "def estimate_relative_humidity(q,e_s,p):\n",
    "    '''estimates relative humidity using hypsometric equation for pressure, virtual temperature and avg temperature'''\n",
    "    import numpy as np\n",
    "    w = q ## approximating the mixing ratio as the specific humidity\n",
    "    w_s = 0.622 * (e_s / p)\n",
    "    RH = 100.0 * (w / w_s)\n",
    "    return (RH)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "IMPORT BETA VALUES AND APPLY EMPIRICAL TRANSFORM "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "direc = '/raid/gergel'\n",
    "pdur_file = 'pduration.nc'\n",
    "pdur_full = xray.open_dataset(os.path.join(direc,pdur_file)) ## pdur beta parameter for John's transform from Matt Jolly\n",
    "## select out gridcell for pdur b parameter\n",
    "lon_ind_1 = np.argmin(np.abs(pdur_full.lon - lon))\n",
    "lat_ind_1 = np.argmin(np.abs(pdur_full.lat - lat))\n",
    "beta = pdur_full.isel_points(lon=[lon_ind_1],lat=[lat_ind_1])\n",
    "\n",
    "pptamt_in = pptamt * 0.04 ## convert from mm to inches \n",
    "pptdur_calc_df = pptamt_in['precipitation_amount'].to_series().map(lambda x: np.round(24 * (1 - (np.exp(-beta['pdur']*x)))) \n",
    "                                                            if (x > 0) else 0) \n",
    "pptdur = xray.DataArray.from_series(pptdur_calc_df).to_dataset()['precipitation_amount'].values[0,:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ITERATE AND CALCULATE 100-HR AND 1000-HR DFM WITH DERIVED RELATIVE HUMIDITY (FROM SPECIFIC HUMIDITY) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## get elevation of grid cell\n",
    "soil_file = '/raid9/gergel/agg_snowpack/soil_avail.txt'\n",
    "elev_corr_info = import_gridcell_elevation(soil_file) \n",
    "h = get_elev_for_lat_lon(elev_corr_info,lat,lon)\n",
    "\n",
    "## get pressure \n",
    "p = estimate_p(h)\n",
    "\n",
    "tmois=np.zeros(shape=(x,7))\n",
    "bv=np.zeros(shape=(x,7))\n",
    "ymc=np.zeros(shape=(x,1))\n",
    "ndays = len(rhmax)\n",
    "\n",
    "fm1000_rh = np.ndarray(shape=(x,ndays),dtype='float')\n",
    "fm100_rh = np.ndarray(shape=(x,ndays),dtype='float')\n",
    "\n",
    "for day in np.arange(ndays):\n",
    "    t_avg = (tmax[day] + tmin[day])/2.0 ## in Kelvin still \n",
    "    e_s = estimate__e_s(t_avg) ## sat vapor pressure\n",
    "    satvpx = estimate__e_s(tmax[day])\n",
    "    satvpn = estimate__e_s(tmin[day])\n",
    "    RH = estimate_relative_humidity(q[day],e_s,p)\n",
    "    ambvp = (RH * e_s) / 100.0 \n",
    "    rhmax = 100.0 * (ambvp/satvpn)\n",
    "    ## constrain RH to be 100% or less \n",
    "    if rhmax > 100:\n",
    "        rhmax = 100\n",
    "    rhmin = 100.0 * (ambvp/satvpx)\n",
    "    if rhmin > 100:\n",
    "        rhmin = 100\n",
    "    tmois,fm1000_rh[0,day],fm100_rh[0,day],bv = calc_fm100_fm1000(pptdur[day],rhmax,rhmin,\n",
    "                                                            kelvin_to_fahrenheit(tmax[day]),\n",
    "                                                            kelvin_to_fahrenheit(tmin[day]),\n",
    "                                                            lat,tmois,bv,julians[day],ymc)\n",
    "    ymc=fm100[0,day]\n",
    "    print(day)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "CREATE ARRAY FOR DATAFRAME INDEX "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "if (scenario == \"historical\"):\n",
    "    tp_start = dt.datetime(1950,1,1)\n",
    "    tp_end = dt.datetime(2006,1,1)\n",
    "else:\n",
    "    tp_start = dt.datetime(2006,1,1)\n",
    "    tp_end = dt.datetime(2100,1,1)\n",
    "\n",
    "    dates = [tp_start + dt.timedelta(days=i) for i in range(0, (tp_end - tp_start).days)]\n",
    "dates_arr = np.asarray(dates)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "WRITE TO NETCDF FILE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "d = {'fm100':fm100[0,:],'fm1000':fm1000[0,:]}\n",
    "df = pd.DataFrame(data=d,index=dates_arr)\n",
    "ds = xray.Dataset.from_dataframe(df)\n",
    "direc = '/raid/gergel/dfm_results/%s_%s' % (model,scenario)\n",
    "if not os.path.exists(direc):\n",
    "    os.makedirs(direc_gc) ## if directory doesn't exist, create it \n",
    "## save to netcdf \n",
    "filename = 'dfm_%s_%s' % (model,scenario,lat,lon)\n",
    "ds.to_netcdf(os.path.join(direc,filename))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "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.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}