_init__.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numpy.ma as ma\n",
    "import pandas as pd\n",
    "import sys\n",
    "import os  \n",
    "import datetime\n",
    "from netCDF4 import Dataset  \n",
    "from osgeo import gdal, ogr\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import NC_CLIP_VIA_SHP\n",
    "import NSE_Func\n",
    "import Get_NMME_index_order\n",
    "import Get_NMME_Var_Dict\n",
    "import Get_NMME_Precip\n",
    "import Bayesian_weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "NCInputBase = ['http://iridl.ldeo.columbia.edu/SOURCES/.Models/.NMME/.CanCM4i/.FORECAST/.MONTHLY/.prec/dods',\n",
    "               'http://iridl.ldeo.columbia.edu/SOURCES/.Models/.NMME/.Cansips/.FORECAST/.MONTHLY/.prec/dods',\n",
    "               'http://iridl.ldeo.columbia.edu/SOURCES/.Models/.NMME/.CanSIPSv2/.FORECAST/.MONTHLY/.prec/dods',\n",
    "               'http://iridl.ldeo.columbia.edu/SOURCES/.Models/.NMME/.CMC1-CanCM3/.FORECAST/.MONTHLY/.prec/dods',\n",
    "               'http://iridl.ldeo.columbia.edu/SOURCES/.Models/.NMME/.CMC2-CanCM4/.FORECAST/.MONTHLY/.prec/dods',\n",
    "               'http://iridl.ldeo.columbia.edu/SOURCES/.Models/.NMME/.COLA-RSMAS-CCSM3/.MONTHLY/.prec/dods',\n",
    "               'http://iridl.ldeo.columbia.edu/SOURCES/.Models/.NMME/.COLA-RSMAS-CCSM4/.MONTHLY/.prec/dods',\n",
    "               'http://iridl.ldeo.columbia.edu/SOURCES/.Models/.NMME/.CPC-CMAP/.prate/dods',\n",
    "               'http://iridl.ldeo.columbia.edu/SOURCES/.Models/.NMME/.CPC-CMAP-URD/.prate/dods',\n",
    "               'http://iridl.ldeo.columbia.edu/SOURCES/.Models/.NMME/.CPC-PRECIP/.prate/dods',\n",
    "               'http://iridl.ldeo.columbia.edu/SOURCES/.Models/.NMME/.GEM-NEMO/.FORECAST/.MONTHLY/.prec/dods',\n",
    "               'http://iridl.ldeo.columbia.edu/SOURCES/.Models/.NMME/.GFDL-CM2p1/.MONTHLY/.prec/dods',\n",
    "               'http://iridl.ldeo.columbia.edu/SOURCES/.Models/.NMME/.GFDL-CM2p1-aer04/.MONTHLY/prec/dods',\n",
    "               'http://iridl.ldeo.columbia.edu/SOURCES/.Models/.NMME/.GFDL-CM2p5-FLOR-A06/.MONTHLY/.prec/dods',\n",
    "               'http://iridl.ldeo.columbia.edu/SOURCES/.Models/.NMME/.GFDL-CM2p5-FLOR-B01/.MONTHLY/.prec/dods',\n",
    "               'http://iridl.ldeo.columbia.edu/SOURCES/.Models/.NMME/.IRI-ECHAM4p5-AnomalyCoupled/.MONTHLY/.prec/dods',\n",
    "               'http://iridl.ldeo.columbia.edu/SOURCES/.Models/.NMME/.IRI-ECHAM4p5-DirectCoupled/.MONTHLY/.prec/dods',\n",
    "               'http://iridl.ldeo.columbia.edu/SOURCES/.Models/.NMME/.NASA-GEOSS2S/.FORECAST/.MONTHLY/.prec/dods',\n",
    "               'http://iridl.ldeo.columbia.edu/SOURCES/.Models/.NMME/.NASA-GMAO/.MONTHLY/.prec/dods',\n",
    "               'http://iridl.ldeo.columbia.edu/SOURCES/.Models/.NMME/.NASA-GMAO-062012/.MONTHLY/.prec/dods',\n",
    "               'http://iridl.ldeo.columbia.edu/SOURCES/.Models/.NMME/.NCAR-CESM1/.FORECAST/.MONTHLY/.prec/dods',\n",
    "               'http://iridl.ldeo.columbia.edu/SOURCES/.Models/.NMME/.NCEP-CFSv1/.MONTHLY/.prec/dods',\n",
    "               'http://iridl.ldeo.columbia.edu/SOURCES/.Models/.NMME/.NCEP-CFSv2/.FORECAST/.EARLY_MONTH_SAMPLES/.MONTHLY/.prec/dods']\n",
    "ShpInputBase = ['C:/Users/Lujun/Desktop/ShapeFIle/CITY_SHP/CHICAGO_IL.shp',\n",
    "                'C:/Users/Lujun/Desktop/ShapeFIle/CITY_SHP/DALLAS_TX.shp',\n",
    "                'C:/Users/Lujun/Desktop/ShapeFIle/CITY_SHP/DENVER_CO.shp',\n",
    "                'C:/Users/Lujun/Desktop/ShapeFIle/CITY_SHP/SIOUXFALL_SD.shp',\n",
    "                'C:/Users/Lujun/Desktop/ShapeFIle/CITY_SHP/OKC_OK.shp',\n",
    "                'C:/Users/Lujun/Desktop/ShapeFIle/CITY_SHP/LA_CA.shp',\n",
    "                'C:/Users/Lujun/Desktop/ShapeFIle/CITY_SHP/ORLA_FL.shp',\n",
    "                'C:/Users/Lujun/Desktop/ShapeFIle/CITY_SHP/PHIL_PA.shp',\n",
    "                'C:/Users/Lujun/Desktop/ShapeFIle/CITY_SHP/SEA_WA.shp']\n",
    "NMME_Member_Name = ['CanCM4i','Cansips','CanSIPSv2','CMC1-CanCM3',\n",
    "                    'CMC1-CanCM4','COLA-RMSAS-CCSM3','COLA-RMSAS-CCSM4',\n",
    "                    'CPC-CMAP','CPC-CMAP-URD','CPC-PRECIP','GEM-NEMO',\n",
    "                    'GFDL-CM2p1','GDFL-CM2p1-aer04','GFDL-CM2P5-FLOR-A06',\n",
    "                    'GFDL-CM2P5-FLOR-B01','IRI-ECHAM4p5-AnomalyCoupled',\n",
    "                    'IIRI-ECHAM4p5-DirectCoupled','NASA-GEOSS2S','NASA-GMAO',\n",
    "                    'NASA-GMAO-062012','NCAR-CESM1','NCEP-CFSv1','NCEP-CFSv2']\n",
    "City_Name = ['Chicago','Dallas','Denver','Sioux_Fall','Oklahoma_City','Los_Angeles','Orlando','Philadaphia','Seattle']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(np.size(NCInputBase)):\n",
    "    for j in range(np.size(ShpInputBase)):\n",
    "        print(NMME_Member_Name[i],City_Name[j])\n",
    "        Out_Direc = '../Proj_NMME/NMMEClipped/'+str(i)+'_'+NMME_Member_Name[i]+'_'+City_Name[j]+'.nc'\n",
    "        ncInput = Dataset(NCInputBase[i])\n",
    "        NMME_Precip_Clip_Via_Shp(ncInput,ShpInputBase[j],Out_Direc,0)\n",
    "        ncInput.close()   # close the input file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Simple averaging for each member of NMME at each locations\n",
    "for i in range(np.size(NMME_Member_Name)):\n",
    "    for j in range(np.size(City_Name)):\n",
    "        print(City_Name[j],NMME_Member_Name[i])\n",
    "        Out_Direc = 'NMMEClipped/'+str(i)+'_'+NMME_Member_Name[i]+'_'+City_Name[j]+'_SMA.nc'\n",
    "        Direc = 'NMMEClipped/'+str(i)+'_'+NMME_Member_Name[i]+'_'+City_Name[j]+'.nc'\n",
    "        ncInput = Dataset(Direc)\n",
    "        var_list = list(ncInput.variables.keys())\n",
    "        if np.size(var_list) == 5:\n",
    "            data = ncInput['prec'][:]\n",
    "            reftime = ncInput['S'][:]\n",
    "            Long = ncInput['X'][:]\n",
    "            Lat = ncInput['Y'][:]\n",
    "            data_tab = data.mean(axis=1)\n",
    "        elif np.size(var_list) == 4:\n",
    "            data = ncInput['prate'][:]\n",
    "            reftime = ncInput['S'][:]\n",
    "            Long = ncInput['X'][:]\n",
    "            Lat = ncInput['Y'][:] \n",
    "        ncOutput = Dataset(Out_Direc, 'w', format='NETCDF4')\n",
    "        ncOutput.createDimension('X', np.size(Long))\n",
    "        ncOutput.createDimension('Y', np.size(Lat))\n",
    "        ncOutput.createDimension('S', np.size(reftime))\n",
    "        # Add lat Variable\n",
    "        var_out_Y = ncOutput.createVariable('Y','f',(\"Y\"))\n",
    "        ncOutput.variables['Y'][:] = Lat[:]\n",
    "        # Add lon Variable\n",
    "        var_out_X = ncOutput.createVariable('X','f',(\"X\"))\n",
    "        ncOutput.variables['X'][:] = Long[:]   \n",
    "        # Add reftime Variable\n",
    "        var_out_S = ncOutput.createVariable('S','f',(\"S\"))\n",
    "        ncOutput.variables['S'][:] = reftime[:]\n",
    "        # Add data Variable\n",
    "        if np.size(var_list) == 5:\n",
    "            var_out_data = ncOutput.createVariable('prec', 'f',(\"S\",\"Y\",\"X\"))\n",
    "            for k in range(np.size(reftime)):\n",
    "                ncOutput.variables['prec'][k,:,:] = data_tab[k,:,:]\n",
    "        else:\n",
    "            var_out_data = ncOutput.createVariable('prate', 'f',(\"S\",\"Y\",\"X\"))\n",
    "            for k in range(np.size(reftime)):\n",
    "                ncOutput.variables['prate'][k,:,:] = data[k,:,:]\n",
    "        print(np.shape(data_tab))\n",
    "        # attr\n",
    "        ncOutput.history = \"CLIP Created datatime \" + datetime.datetime.now().strftime(\"%Y-%m-%d %H:%M:%S\") + \" by LujunZ at OU\"\n",
    "        ncOutput.source  = \"netCDF4 under python 3.6.5\"\n",
    "        ncOutput.close()  # close the new file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read and Extract reference data\n",
    "RefNetCDFInput = 'C:/Users/Lujun/Desktop/2019 Fall Semester/20191015 NMME_Project/precip.mon.total.1x1.v2018.nc'\n",
    "ncInput = Dataset(RefNetCDFInput)\n",
    "for i in range(len(City_Name)):\n",
    "    Out_Direc = 'RefClipped/'+City_Name[i]+'_Ref_Data.nc'\n",
    "    shpDS = ogr.Open(ShpInputBase[i])\n",
    "    shpLyr = shpDS.GetLayer()\n",
    "    Envelop = shpLyr.GetExtent() \n",
    "    xmin,ymin,xmax,ymax = [int(Envelop[0]),\n",
    "                           int(Envelop[2]),\n",
    "                           int(Envelop[1]),\n",
    "                           int(Envelop[3])] #Your extents as given above\n",
    "    lon_Ori = ncInput['lon'][:]-180\n",
    "    lat_Ori = ncInput['lat'][:]\n",
    "    # Reference diff\n",
    "    # lat_Ori = lat_Ori\n",
    "    reftime_Ori = ncInput['time'][:]\n",
    "    refData_Ori = ncInput['precip'][:]\n",
    "    #get boundary and xs ys\n",
    "    lat_bnds, lon_bnds = [ymin, ymax], [xmin, xmax]\n",
    "    lat_inds = np.where((lat_Ori >= lat_bnds[0]) & (lat_Ori <= lat_bnds[1]))\n",
    "    lon_inds = np.where((lon_Ori >= lon_bnds[0]) & (lon_Ori <= lon_bnds[1]))\n",
    "    ncols = len(lon_inds[0])\n",
    "    nrows = len(lat_inds[0])\n",
    "    nreftime = len(reftime_Ori)\n",
    "    #create geotransform\n",
    "    if ncols>0:\n",
    "        xres = (xmax - xmin) / float(ncols)\n",
    "    else:\n",
    "        xres = 1 \n",
    "    if nrows>0:\n",
    "        yres = (ymax - ymin) / float(nrows)\n",
    "    else:\n",
    "        yres = 1\n",
    "    geotransform = (xmin,xres,0,ymax,0, -yres)\n",
    "    #create mask\n",
    "    #mask_DS = gdal.GetDriverByName('MEM').Create('', ncols, nrows, 1 ,gdal.GDT_Int32)\n",
    "    #mask_RB = mask_DS.GetRasterBand(1)\n",
    "    #mask_RB.Fill(0) #initialise raster with zeros\n",
    "    #mask_RB.SetNoDataValue(-32767)\n",
    "    #mask_DS.SetGeoTransform(geotransform)\n",
    "    #maskvalue = 1\n",
    "    #err = gdal.RasterizeLayer(mask_DS, [maskvalue], shpLyr)\n",
    "    #mask_DS.FlushCache()\n",
    "    #mask_array = mask_DS.GetRasterBand(1).ReadAsArray()    \n",
    "    #mask_RES = ma.masked_equal(mask_array, 255)          \n",
    "    #ma.set_fill_value(mask_RES, -32767)  \n",
    "    # Subset\n",
    "    ref_var_subset = refData_Ori[:,min(lat_inds[0])-1:max(lat_inds[0]),min(lon_inds[0])-1:max(lon_inds[0])]\n",
    "    #if mask_RES.mask:\n",
    "    #    ref_var_subset.__setmask__(np.logical_not(np.flipud(mask_RES.mask)))  # update mask (flipud is reverse 180)\n",
    "    lon_subset = lon_Ori[lon_inds]\n",
    "    lat_subset = lat_Ori[lat_inds]\n",
    "    \n",
    "    ncOutput = Dataset(Out_Direc, 'w', format='NETCDF4')\n",
    "\n",
    "    #ncOutput.createDimension('time', None)\n",
    "    ncOutput.createDimension('X', ncols)\n",
    "    ncOutput.createDimension('Y', nrows)\n",
    "    ncOutput.createDimension('S', nreftime)\n",
    "\n",
    "    # Add lat Variable\n",
    "    var_out_Y = ncOutput.createVariable('Y','f',(\"Y\"))\n",
    "    ncOutput.variables['Y'][:] = lat_subset[:]\n",
    "\n",
    "    # Add lon Variable\n",
    "    var_out_X = ncOutput.createVariable('X','f',(\"X\"))\n",
    "    ncOutput.variables['X'][:] = lon_subset[:]\n",
    "\n",
    "    # Add reftime Variable\n",
    "    var_out_S = ncOutput.createVariable('S','f',(\"S\"))\n",
    "    ncOutput.variables['S'][:] = reftime_Ori[:]\n",
    "\n",
    "    # Add data Variable\n",
    "    var_out_data = ncOutput.createVariable('prec', 'f',(\"S\",\"Y\",\"X\"))\n",
    "    for k in range(nreftime):\n",
    "        ncOutput.variables['prec'][k,:,:] =  ref_var_subset[k,:,:]\n",
    "    #print(sum(ref_var_subset))\n",
    "    # attr\n",
    "    ncOutput.history = \"CLIP Created datatime\" + datetime.datetime.now().strftime(\"%Y-%m-%d %H:%M:%S\") + \" by LujunZ at OU\"\n",
    "    ncOutput.source  = \"netCDF4 under python 3.6.5\"\n",
    "    ######################################################\n",
    "    #                   Write close                      #\n",
    "    ######################################################\n",
    "    ncOutput.close()  # close the new file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Simple averaging all members of NMME at each location\n",
    "# Create empty Pandas Dataframe \n",
    "# Montly Data from Jan 1st 1980 to Dec 31st 2016\n",
    "\n",
    "# Overall_Date_Range = pd.date_range(start=\"18910101\", end=\"20191201\", freq=\"M\")\n",
    "NMME_Date_Range = pd.date_range(start=\"19600101\", end=\"20191130\", freq=\"M\") \n",
    "Ref_Date_Range = pd.date_range(start=\"18910101\", end=\"20161231\", freq=\"M\") #Montly Data from 1982 to 2010\n",
    "# Change data range here\n",
    "Selected_Date_Range = pd.date_range(start=\"19800101\", end=\"20161231\", freq=\"M\") \n",
    "NMME_Data = pd.DataFrame(index =NMME_Date_Range, columns=NMME_Member_Name)\n",
    "\n",
    "# Read all data \n",
    "for i in range(np.size(City_Name)):\n",
    "    print(City_Name[i])\n",
    "    NMME_Data = pd.DataFrame(index =Selected_Date_Range, columns=NMME_Member_Name)\n",
    "    Out_Direc = 'SMA/'+City_Name[i]+'_SMA.xlsx'\n",
    "    ncRefInput = Dataset('RefClipped/'+City_Name[i]+'_Ref_Data.nc')\n",
    "    refData = ncRefInput['prec'][:]\n",
    "    #print(sum(refData))\n",
    "    refLength = np.size(ncRefInput['S'][:])\n",
    "    tab_ref_data = np.zeros(refLength)\n",
    "    for m in range(refLength):\n",
    "        tab_ref_data[m]=refData[m,:,:].mean()/30  \n",
    "    tab_ref_data = pd.Series(tab_ref_data,index=Ref_Date_Range)\n",
    "    tab_ref_data = tab_ref_data[\"19800101\":\"20161231\"]\n",
    "    \n",
    "    for j in range(np.size(NMME_Member_Name)):\n",
    "        Direc = 'NMMEClipped/'+str(j)+'_'+NMME_Member_Name[j]+'_'+City_Name[i]+'_SMA.nc'\n",
    "        ncInput = Dataset(Direc)\n",
    "        Var_name = list(ncInput.variables.keys())\n",
    "        if 'prec' in Var_name:\n",
    "            data = ncInput['prec'][:]\n",
    "        else:\n",
    "            data = ncInput['prate'][:]\n",
    "        reftime = ncInput['S'][:]\n",
    "        tab_idx = pd.date_range(NMME_Date_Range[int(reftime[0])], \n",
    "                                periods=np.size(reftime), freq='M')        \n",
    "        tab_data = np.zeros(len(reftime))\n",
    "        for k in range(len(reftime)):\n",
    "            tab_data[k]=data[k,:,:].mean()\n",
    "            tab_data = pd.Series(tab_data,index=tab_idx)\n",
    "        NMME_Data[NMME_Member_Name[j]]=tab_data\n",
    "        NMME_Data['SMA'] = NMME_Data.mean(1)\n",
    "        NMME_Data['RefData'] = tab_ref_data\n",
    "    NMME_Data.to_excel(Out_Direc, sheet_name='Sheet_name_1')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Data reading from former excel \n",
    "# Data cleaning using only selected member\n",
    "drop_name = ['CanCM4i','Cansips','CanSIPSv2','CMC1-CanCM3','CMC1-CanCM4','GEM-NEMO',\n",
    "             'GFDL-CM2P5-FLOR-A06','NASA-GEOSS2S','NASA-GMAO','NASA-GMAO-062012',\n",
    "             'NCAR-CESM1','NCEP-CFSv1','NCEP-CFSv2']\n",
    "for i in range(np.size(City_Name)):\n",
    "    Xlsx_name = 'SMA/'+City_Name[i]+'_SMA.xlsx'\n",
    "    exec('data_'+City_Name[i]+'=pd.read_excel(\"'+Xlsx_name+'\", sheet_name=\"Sheet_name_1\",index_col=0)')\n",
    "    for j in range(np.size(drop_name)):\n",
    "        exec('data_'+City_Name[i]+'=data_'+City_Name[i]+'.drop([\"'+drop_name[j]+'\"],axis=1)')\n",
    "    # here specify the study period\n",
    "    #exec('data_'+City_Name[i]+'=data_'+City_Name[i])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Computing the Bayesian weights for different members at each locations\n",
    "# Calculating the Bayesian Model Averaging at each loceations\n",
    "# Saving the corresponding results\n",
    "for i in range(np.size(City_Name)):\n",
    "    Out_Direc = 'BMA/'+City_Name[i]+'_BMA.xlsx'\n",
    "    exec('data = data_'+City_Name[i]+'.drop([\"RefData\"],axis=1)')\n",
    "    exec('data_'+City_Name[i]+'[\"SMA\"]=data.mean(1)')\n",
    "    exec('ref = data_'+City_Name[i]+'[\"RefData\"]')\n",
    "    ref_inp = np.array(ref[\"19820101\":\"19991231\"])\n",
    "    data_inp = np.array(data[\"19820101\":\"19991231\"])\n",
    "    exec('weights_'+City_Name[i]+'=Bayesian_weights(data_inp,ref_inp)')\n",
    "    exec('data_'+City_Name[i]+'[\"BMA\"] = data_'+City_Name[i]+'\\\n",
    "         .drop([\"RefData\",\"SMA\"],axis=1).values.dot(np.reshape(weights_'+City_Name[i]+',-1,1))')\n",
    "    exec('data_'+City_Name[i]+'.to_excel(Out_Direc, sheet_name=\"Sheet_name_1\")') "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}