%matplotlib inline
import pandas as pd
import zipfile
import shutil
import urllib2

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.dates as dates


from urllib2 import urlopen     
import time

from bs4 import BeautifulSoup
from pylab import rcParams
import platform
rcParams['figure.figsize'] = 15, 10
import re
import os

import glob
import urllib

import gzip

import ftplib
import calendar
import datetime
from datetime import date

import pymodis

import arcpy
arcpy.CheckOutExtension("spatial")
from arcpy.sa import *

print("Operating System " + platform.system() + " " + platform.release())
print("Python Version " + str(sys.version))
print("Pandas Version " + str(pd.__version__))
print("Numpy Version " + str(np.__version__))

import UBM
UBM.__version__

"""
Explanation: Import Libraries
End of explanation
"""


#Download HDF Files
#'h09v05','h09v04','h08v05'
tiles = ['h09v05','h09v05','h09v04','h08v05']
save_path = 'H:/GIS/MODIS/'

UBM.get_modis(tiles, save_path)

files = UBM.get_file_list(save_path)

"""
Explanation: MODIS16
Potential and Actual Evapotranspiration
https://github.com/lucadelu/pyModis/blob/master/docs/source/examples/pyModis.ipynb<br>
http://www.ntsg.umt.edu/project/mod16
Download HDF Files
End of explanation
"""


save_path = 'H:/GIS/MODIS/'
UBM.reproject_modis(files,save_path,'ET')

path = "H:/GIS/MODIS/ET/"
outpath="H:/GIS/MODIS16.gdb/"
UBM.clip_and_fix(path,outpath,'ET')

"""
Explanation: Reproject MODIS Data
These scripts reproject the hdf files from the <a href="https://modis-land.gsfc.nasa.gov/MODLAND_grid.html">wacky sinusoidal MODIS projection</a> to <a href="http://spatialreference.org/ref/epsg/4269/">NAD83 Zone 12</a> for analysis.  In this section, the MODIS rasters are also clipped to Utah Watersheds (see image below) and the fill values are made to null values. Fill values are described in the <a href="http://files.ntsg.umt.edu/data/NTSG_Products/MOD16/MOD16_global_evapotranspiration_description.pdf">MODIS16 documentation</a>:
<br>
<ul>
<li>Fill value, out of the earth 32767</li>
<li>Water body 32766</li>
<li>Barren or sparsely vegetated 32765</li> 
<li>Permanent snow and ice 32764</li>
<li>Permanent wetland 32763</li>
<li>Urban or Built-up 32762</li>
<li>Unclassified 32761</li>
<img src="https://cfpub.epa.gov/surf/images/states/ut.gif">

### Reproject ET
End of explanation
"""


save_path = "H:/GIS/MODIS/PET/"
UBM.reproject_modis(files,save_path,'PET')

path = "H:/GIS/MODIS/PET/"
outpath="H:/GIS/MODIS16.gdb/"
UBM.clip_and_fix(path,outpath,'PET')

"""
Explanation: Reproject PET
End of explanation
"""


path="H:/GIS/MODIS16.gdb/"
UBM.merge_rasts(path,data_type='ET',monthRange=[1,12],yearRange=[2000,2015])

"""
Explanation: Merge Data
MODIS data is downloaded as separate cells based on the Sinusoidal grid.  Three MODIS cells cover Utah ('h09v05','h09v04','h08v05'). The following scripts mosaic (merge) the three individual rasters for each month into one seamless monthly raster for the entire state. 
<img src="https://modis-land.gsfc.nasa.gov/images/MODIS_sinusoidal_grid1.gif" alt="MODIS grid">
Merge ET
End of explanation
"""


path="H:/GIS/MODIS16.gdb/"
UBM.merge_rasts(path,data_type='PET',monthRange=[1,12],yearRange=[2000,2015])

"""
Explanation: Merge PET
End of explanation
"""


path="H:/GIS/MODIS16.gdb/"
arcpy.env.workspace = path

scalingFactor = 10000.0 #convert to m/month
out_path="H:/GIS/MODIS16.gdb/"
UBM.scale_modis(path, out_path, scaleby = 10000.0, data_type = 'ET')
UBM.scale_modis(path, out_path, scaleby = 10000.0, data_type = 'PET')

"""
Explanation: Scale MODIS Data
The documentation for the dataset says that the raw files have to be multiplied by 0.1 to get mm/month.  To convert to meters per month, we have to multiply the raw files by 0.0001 (0.1 x 0.001) or divide by 10,000.  These scripts do that division.
End of explanation
"""


path="H:/GIS/MODIS16/"
arcpy.env.workspace = path

for rast in arcpy.ListRasters('*c'):
    print(rast)
    arcpy.Delete_management(rast, 'RasterDataset')
    
for rast in arcpy.ListRasters('*h*v*'):
    print(rast)
    arcpy.Delete_management(rast, 'RasterDataset')

"""
Explanation: Clean up GDB
This deletes intermediate files left over from previous processing steps. For whatever reason, it take a while.
End of explanation
"""


path= "H:/GIS/MODIS2/MODIS.gdb/"
outpath = "H:/GIS/MODIS2/MODIS16.gdb/"

units = 'CELL'
radius = 15

arcpy.env.workspace = path

def fill_holes(path, outpath, wildcard, units='CELL', radius=15):
    for rast in arcpy.ListRasters(wildcard):
        rfilled = arcpy.sa.Con(arcpy.sa.IsNull(rast),
                              arcpy.sa.FocalStatistics(rast,
                                                       arcpy.sa.NbrCircle(radius, units),'MEAN'), rast)
        dsc = arcpy.Describe(rast)
        nm = dsc.baseName
        rfilled.save(outpath+nm)
        print(nm)
    

"""
Explanation: Fill Holes in Rasters
This processing step fills in null values created when the fill values were removed above.
From: http://gis.stackexchange.com/questions/136075/fill-in-nodata-gaps-in-raster-using-arcgis-for-desktop
End of explanation
"""


L = u'H:\GIS\SNODAS'
UBM.get_snodas(L)

"""
Explanation: SNODAS
Data were downloaded using polaris: http://nsidc.org/data/polaris/
SNODAS ftp Download
Alternatively, you can download from the ftp, which is slower. <br>
ftp://sidads.colorado.edu/pub/DATASETS/NOAA/G02158/ <br>
https://support.nsidc.org/entries/64231694-FTP-Client-Data-Access 
Download SNODAS TARs
End of explanation
"""


#https://github.com/PSU-CSAR/SNODAS-SWE/blob/master/snodas-swe-prep.py

"""
Explanation: Extract SNODAS Data Files
http://nsidc.org/data/docs/noaa/g02158_snodas_snow_cover_model/
End of explanation
"""


snodasTars = glob.glob(u'H:/GIS/SNODAS/*.tar')
for tared in snodasTars:
    untar(tared,u'H:/GIS/SNODAS/SNODASUNZIPPED/')

"""
Explanation: Unzip tar files and put in unzip file
End of explanation
"""


unsnodasTars = glob.glob(u'H:/GIS/SNODAS/SNODASUNZIPPED/*.gz')
for tared in unsnodasTars:
    ungz(tared, deletesource=True)

"""
Explanation: Ungz gz files and delete gz files
End of explanation
"""


for hdrfile in glob.glob("H:/GIS/SNODAS/SNODASUNZIPPED/*.Hdr"):
    replace_hdr_file(hdrfile)

"""
Explanation: Replace header file to make compatible with ArcGIS
End of explanation
"""


path = "C:/Users/PAULINKENBRANDT/Downloads/NSIDC_Data"           
UBM.rename_polaris_snodas(path)

"""
Explanation: Polaris Download
Downloaded data from http://nsidc.org/data/polaris/, which allows for clipping to an area, as well as output as GeoTiff Format.  The download takes about 3 days and the unzipping takes about 2 hours.
<b>Extent:</b><br>
<ul type="disc">
<li>Top `44.0254000046`</li>
<li>Left `-116.075366662`</li>
<li>Right `-108.375366657`</li>
<li>Bottom `36.0087333333`</li>

<b>Selected Variables:</b> <p>Snow Melt Runoff at the Base of the Snow Pack, Sublimation of Blowing Snow, Snow Pack Average Temperature, Sublimation from the Snow Pack, Snow Water Equivalent, Liquid Precipitation, Snow Depth, Solid Precipitation</p><br>
<b>Download Format:</b> GeoTIFF<br>
<b>Output Grid:</b> 30" (Cylindrical Equidistant)

After unzipping files, rename them to make them easier to handle.

From http://nsidc.org/data/docs/noaa/g02158_snodas_snow_cover_model/, the file abbreviations are as follows:
<ul type="disc">
<li>`RAIN` = `Wet Precip`</li>
<li>`SWEQ` = `Snow Water Equivalent`</li> 
<li>`SNOD` = `Snow Depth`</li> 
<li>`SPAT` = `Snow Pack Average Temp`</li>
<li>`BSSB` = `Blowing Snow Sublimation`</li>
<li>`SNML` = `Snowmelt`</li>
<li>`SPSB` = `Snow Pack Sublimation`</li>
</ul>
End of explanation
"""


UBM.snow_summary('RAIN',10000.0, monthRange = [1,12], yearRange = [2015,2015])

"""
Explanation: Data Merge (daily to monthly)
This merge uses geoTiff files from Polaris
https://pro.arcgis.com/en/pro-app/help/analysis/spatial-analyst/mapalgebra/building-complex-statements.htm
Precipitation
Includes snow and rain, which are provided as separate data sets in the SNOTEL data.
Wet Precipitation
This is a measure of the cumulative incoming rain over a month.<br> 
Units are meters of water per month.
End of explanation
"""


UBM.snow_summary('SNOW',10000.0)

"""
Explanation: Snow
This represents the "dry" precipitation that falls to the ground as snow.<br>
We summed daily data to create cumulative monthly snowfall.<br>
Units are in meters of water per month.
End of explanation
"""


monthRange = [1,12] 
yearRange = [2003,2016]

g = {}
path="H:/GIS/SNODAS.gdb/"
arcpy.env.workspace = path
arcpy.env.overwriteOutput = True

for y in range(yearRange[0],yearRange[1]+1): #set years converted here
    for m in range(monthRange[0],monthRange[1]+1): #set months converted here
        my = str(y)+str(m).zfill(2)
        newdn = 'TPPT' + my
        try:
            calc = Plus('SNOW'+ my +'SUM', 'RAIN'+ my +'SUM')
            calc.save(newdn+'SUM')
            print(newdn)
        except(RuntimeError):
            pass

"""
Explanation: Total Precipitation
End of explanation
"""


UBM.snow_summary('SWEQ', 1000.0, 'MEDIAN', monthRange=[1, 12], yearRange=[2003, 2016])

"""
Explanation: SWE
SWE = Snow-water equivalent;  This is a measure of the estimated total water content of the snow pack. <br>
The median is used to summarize these data.<br>
Units are in meters of water.
End of explanation
"""


UBM.snow_summary('SNML', 100000.0, 'SUM', monthRange=[1, 12], yearRange=[2003, 2015])

"""
Explanation: Snowmelt
SNML = Cumulative monthly snowmelt calculated by summing daily snowmelt data.<br>
Units are in meters of water per month.
Snow Melt Runoff at the Base of the Snow Pack.
Units are meters of water per month.
End of explanation
"""


UBM.snow_summary('BSSB', 100000.0, 'SUM', monthRange=[9, 9], yearRange=[2016, 2016])

UBM.snow_summary('SPSB', 100000.0, 'SUM', monthRange=[9, 9], yearRange=[2016, 2016])

g = {}
path="H:/GIS/SNODAS.gdb/"
arcpy.env.workspace = path
arcpy.env.overwriteOutput = True

monthRange = [1,12] 
yearRange = [2003,2016]

for y in range(yearRange[0],yearRange[1]+1): #set years converted here
    for m in range(monthRange[0],monthRange[1]+1): #set months converted here
        my = str(y)+str(m).zfill(2)
        newdn = 'TSSB' + my
        try:
            calc = Plus('BSSB'+ my +'SUM', 'SPSB'+ my +'SUM')
            calc.save(newdn+'SUM')
            print(newdn)
        except(RuntimeError):
            pass

"""
Explanation: Snow Sublimation
BSSB = Cumulative monthly blowing snow sublimation.<br>
SPSB = Cumulative monthly snow pack sublimation.<br>
TSSB = Total Cumulative monthly snow sublimation.<br>
Units are in meters of water per month.
End of explanation
"""


UBM.totalavg('TSSB')

UBM.totalavg('SNML', monthRange=[2, 12], yearRange=[2003, 2016])

"""
Explanation: Monthly Averages
End of explanation
"""


soildir = 'H:/GIS/soils/STATSGO/'

keeplist = ['hzname', 'hzdept_r','hzdepb_r','hzthk_r','dbovendry_r','ksat_r','awc_r','wthirdbar_r',
            'wfifteenbar_r','mukey','cokey','chkey']

chorizon = pd.read_excel(soildir+'chorizon.xlsx')

droplist = list(set(list(chorizon.columns)) - set(keeplist))
chorizon.drop(droplist, axis=1, inplace=True)

chorizon.columns = [str(x).strip() for x in list(chorizon.columns)]

chorizon.columns

"""
Explanation: Soil Data
<img src='https://casoilresource.lawr.ucdavis.edu/files/6114/1807/4326/SSURGO-analysis-diagram.preview.png'>
STATSGO
End of explanation
"""


chorizon['hzthk_r'] = chorizon['hzdepb_r'] - chorizon['hzdept_r']
thick = chorizon.groupby(['cokey'])['hzthk_r'].sum().to_dict()
chorizon['comthk'] = chorizon['cokey'].apply(lambda x: thick[x],1)
chorizon['thickness_weight'] = chorizon['hzthk_r']/chorizon['comthk']

"""
Explanation: Get thickness of all measured horizons and calculate weighting factor for each horizon.
End of explanation
"""


chorizon['tw_ksat'] = chorizon[['thickness_weight','ksat_r']].apply(lambda x: x[0]*x[1],1)
chorizon['tw_wilt'] = chorizon[['thickness_weight','wfifteenbar_r']].apply(lambda x: x[0]*x[1],1)
chorizon['tw_field_capacity'] = chorizon[['thickness_weight','wthirdbar_r']].apply(lambda x: x[0]*x[1],1)
chorizon['tw_dense'] = chorizon[['thickness_weight','dbovendry_r']].apply(lambda x: x[0]*x[1],1)

"""
Explanation: Calculated weighting of each important horizon variable.
End of explanation
"""


component_props = chorizon.groupby('cokey').sum()
component_props.drop(['hzdept_r','hzdepb_r','dbovendry_r','ksat_r','awc_r','wthirdbar_r','wfifteenbar_r'],axis=1,inplace=True)

component_props.dropna(subset=['tw_ksat','tw_wilt','tw_field_capacity','tw_dense'],how='all',inplace=True)

component_props.reset_index(inplace=True)

"""
Explanation: Combine horizons into components. Drop na and unneeded fields.
End of explanation
"""


component = pd.read_excel(soildir+'component.xlsx')
droplist = list(set(list(component.columns))-set(['comppct_r','compname','mukey','cokey']))
component.drop(droplist, axis=1, inplace=True)

com_props = pd.merge(component_props, component, on='cokey')

corest = pd.read_excel('H:/GIS/soils/STATSGO/corestrictions.xlsx')
corest.drop([u'reskind', u'reshard', u'resdept_l', u'resdept_h', u'resdepb_l', u'resdepb_r', u'resdepb_h', 
             u'resthk_l', u'resthk_r', u'resthk_h', u'corestrictkey'], inplace=True, axis=1)

comprops = pd.merge(corest,com_props, on='cokey')

"""
Explanation: Read in component table. Keep only keys and component percent field.
End of explanation
"""


cowt = comprops.groupby(['mukey'])['comppct_r'].sum().to_dict()
comprops['comppct_tot'] = comprops['mukey'].apply(lambda x: cowt[x],1)
comprops['comppct_r'] = comprops['comppct_r'] / comprops['comppct_tot']

comprops['thickness_m'] = comprops['comppct_r']*comprops['resdept_r']/100.0 # convert to meters
comprops['ksat'] = comprops['tw_ksat'] * comprops['comppct_r']
comprops['wilt'] = comprops['tw_wilt'] * comprops['comppct_r']
comprops['field_cap'] = comprops['tw_field_capacity'] * comprops['comppct_r']
comprops['dense'] = comprops['tw_dense'] * comprops['comppct_r']

muprops = comprops.groupby(['mukey']).sum()

dropfields = [u'cokey', u'hzthk_r', u'chkey', u'comthk', u'tw_ksat', u'tw_wilt',
              u'tw_field_capacity', u'tw_dense', u'thickness_weight','comppct_tot']
muprops.drop(dropfields,axis=1,inplace=True)

muprops['porosity'] = 1 - (muprops['dense']/2.65) # bulk density of material / density of quartz = percent solids
muprops['max_soil_cap_m'] = muprops['thickness_m']*muprops['porosity'] # meters of water
muprops['field_cap_m'] = muprops['max_soil_cap_m']*muprops['field_cap']
muprops['wilt_m'] = muprops['max_soil_cap_m']*muprops['wilt']

muprops.drop(['resdept_r','comppct_r','dense'],axis=1,inplace=True)

muprops.index = muprops.index.astype('str')

muprops.to_pickle('H:/GIS/soils/STATSGO/muprops.pickle')

muprops = pd.read_pickle('H:/GIS/soils/STATSGO/muprops.pickle')

muprops

path = 'H:/GIS/soils/STATSGO/STATSGO.gdb'
arcpy.env.workspace = path
arcpy.env.overwriteOutput = True

fc = 'H:/GIS/soils/STATSGO/STATSGO.gdb/gsmsoilmu_hucs'
narray = muprops.to_records()

arcpy.da.ExtendTable(fc, "MUKEY", narray, "mukey", append_only=False)

path = 'H:/GIS/soils/STATSGO/STATSGO.gdb'
arcpy.env.workspace = path
arcpy.env.overwriteOutput = True

convfields = ['thickness_m','ksat','wilt_m','porosity','max_soil_cap_m','field_cap_m']

for field in convfields:
    arcpy.PolygonToRaster_conversion(fc, field, field, 'CELL_CENTER','', 100)

"""
Explanation: Combine compenents for each mukey.
End of explanation
"""


path = 'H:/GIS/soils/SSURGO/gssurgo_g_ut.gdb'

arcpy.env.workspace = path
arcpy.env.overwriteOutput = True

def arc_to_df(table, fieldlist):
    fields = arcpy.ListFields(table)
    return pd.DataFrame(arcpy.da.TableToNumPyArray(table,fieldlist,null_value=np.nan))

keeplist = ['hzname', 'hzdept_r','hzdepb_r','hzthk_r','dbovendry_r','ksat_r','awc_r','wthirdbar_r',
            'wfifteenbar_r','cokey','chkey']

chorizon = arc_to_df('chorizon',keeplist)

chorizon['hzthk_r'] = chorizon['hzdepb_r'] - chorizon['hzdept_r']
thick = chorizon.groupby(['cokey'])['hzthk_r'].sum().to_dict()
chorizon['comthk'] = chorizon['cokey'].apply(lambda x: thick[x],1)
chorizon['thickness_weight'] = chorizon['hzthk_r']/chorizon['comthk']

chorizon['tw_ksat'] = chorizon[['thickness_weight','ksat_r']].apply(lambda x: x[0]*x[1],1)
chorizon['tw_wilt'] = chorizon[['thickness_weight','wfifteenbar_r']].apply(lambda x: x[0]*x[1],1)
chorizon['tw_field_capacity'] = chorizon[['thickness_weight','wthirdbar_r']].apply(lambda x: x[0]*x[1],1)
chorizon['tw_dense'] = chorizon[['thickness_weight','dbovendry_r']].apply(lambda x: x[0]*x[1],1)

component_props = chorizon.groupby('cokey').sum()
component_props.drop(['hzdept_r','hzdepb_r','dbovendry_r','ksat_r','awc_r','wthirdbar_r','wfifteenbar_r'],axis=1,inplace=True)

component_props.dropna(subset=['tw_ksat','tw_wilt','tw_field_capacity','tw_dense'],how='all',inplace=True)

component_props.reset_index(inplace=True)

component = arc_to_df('component',['comppct_r','compname','mukey','cokey'])

com_props = pd.merge(component_props, component, on='cokey')

fields = arcpy.ListFields('corestrictions')
fieldlist = [field.name for field in fields]
droplist = [u'reskind', u'reshard', u'resdept_l', u'resdept_h', u'resdepb_l', u'resdepb_r', u'resdepb_h', 
             u'resthk_l', u'resthk_r', u'resthk_h', u'corestrictkey']
keepfields = list(set(fieldlist)-set(droplist))

corest = arc_to_df('corestrictions',keepfields)

comprops = pd.merge(corest,com_props, on='cokey')

cowt = comprops.groupby(['mukey'])['comppct_r'].sum().to_dict()
comprops['comppct_tot'] = comprops['mukey'].apply(lambda x: cowt[x],1)
comprops['comppct_r'] = comprops['comppct_r'] / comprops['comppct_tot']

comprops['thickness_m'] = comprops['comppct_r']*comprops['resdept_r']/100.0 # convert to meters
comprops['ksat'] = comprops['tw_ksat'] * comprops['comppct_r']
comprops['wilt'] = comprops['tw_wilt'] * comprops['comppct_r']
comprops['field_cap'] = comprops['tw_field_capacity'] * comprops['comppct_r']
comprops['dense'] = comprops['tw_dense'] * comprops['comppct_r']

muprops = comprops.groupby(['mukey']).sum()

dropfields = [u'hzthk_r', u'comthk', u'tw_ksat', u'tw_wilt',
              u'tw_field_capacity', u'tw_dense', u'thickness_weight','comppct_tot']
muprops.drop(dropfields,axis=1,inplace=True)

muprops['porosity'] = 1 - (muprops['dense']/2.65) # bulk density of material / density of quartz = percent solids
muprops['max_soil_cap_m'] = muprops['thickness_m']*muprops['porosity'] # meters of water
muprops['field_cap_m'] = muprops['max_soil_cap_m']*muprops['field_cap']
muprops['wilt_m'] = muprops['max_soil_cap_m']*muprops['wilt']

muprops.drop(['resdept_r','comppct_r','dense'],axis=1,inplace=True)

muprops.index = muprops.index.astype('str')

muprops.to_pickle('H:/GIS/soils/SSURGO/muprops.pickle')

muprops = pd.read_pickle('H:/GIS/soils/SSURGO/muprops.pickle')

muprops.drop(['OBJECTID'],axis=1,inplace=True)

fc = 'H:/GIS/soils/SSURGO/gssurgo_g_ut.gdb/gmupoly'
narray = muprops.to_records()

fieldlist = list(set(muprops.columns)-set(['mukey']))
arcpy.da.ExtendTable(fc, "MUKEY", narray, "mukey", append_only=False)

path = 'H:/GIS/soils/SSURGO/ssurgo.gdb'

arcpy.env.workspace = path
arcpy.env.overwriteOutput = True

fc = 'gmupoly'
convfields = ['thickness_m','ksat','wilt_m','porosity','max_soil_cap_m','field_cap_m']

for field in convfields:
    arcpy.PolygonToRaster_conversion(fc, field, field, 'CELL_CENTER','', 100)

"""
Explanation: SSURGO
End of explanation
"""


convfields = ['thickness_m','ksat','wilt_m','porosity','max_soil_cap_m','field_cap_m']
path1 = 'H:/GIS/soils/SSURGO/ssurgo.gdb/'
path2 = 'H:/GIS/soils/STATSGO/STATSGO.gdb/'
for field in convfields:
    rfilled = arcpy.sa.Con(arcpy.sa.IsNull(path1+field), path2+field, path1+field)
    rfilled.save('H:/GIS/Soils.gdb/'+field)

"""
Explanation: Combining Soils Data
End of explanation
"""


monthRange = [1,12]

fileplace = 'H:/GIS/PRISM/PRISM'

for m in range(monthRange[0],monthRange[1]+1): 
    testfile = urllib.URLopener()
    testfile.retrieve("http://services.nacse.org/prism/data/public/normals/800m/ppt/"+ str(m).zfill(2), 
                      fileplace+str(m).zfill(2)+'.zip')


testfile = urllib.URLopener()
testfile.retrieve("http://services.nacse.org/prism/data/public/normals/800m/ppt/"+ '14', 
                    fileplace+'14'+'.zip')

fileplace = 'H:/GIS/PRISM/PRISM/'

def unzipper(filepath):

    import zipfile

    f = zipfile.ZipFile(filepath,'r')
    f.extractall(filepath[:-6])
 

for files in glob.glob(fileplace+'*.zip'):
    unzipper(files)
    print(files)

monthRange = [1,12]

fileplace = 'H:/GIS/PRISM/PRISM'

for m in range(monthRange[0],monthRange[1]+1): 
    testfile = urllib.URLopener()
    testfile.retrieve("http://services.nacse.org/prism/data/public/normals/800m/ppt/"+ str(m).zfill(2), 
                      fileplace+str(m).zfill(2)+'.zip')
testfile = urllib.URLopener()
testfile.retrieve("http://services.nacse.org/prism/data/public/normals/800m/ppt/"+ 14, 
                    fileplace+str(m).zfill(2)+'.zip')

"""
Explanation: PRISM
PRISM Data
http://www.prism.oregonstate.edu/
<a href=http://www.prism.oregonstate.edu/documents/PRISM_downloads_web_service.pdf>PRISM Web Services</a>
End of explanation
"""


path="H:/GIS/SNODAS.gdb/"
path2 = "U:/GWP/Groundwater/Projects/BCM/Data/AvailableWater2.gdb/"

def calc_avail_water(path, path2, months='',years='')

    arcpy.env.workspace = path
    arcpy.env.overwriteOutput = True

    if months == '':
        months = [1,12] 
    if years == '':
        years = [2004,2014]

    for y in range(years[0], years[1]+1): #set years converted here
        for m in range(months[0], months[1]+1): #set months converted here
            my = str(y)+str(m).zfill(2)
            newdn = 'AVWT' + my
            rain = Raster('RAIN'+ my +'SUM')
            melt = Raster('SNML'+ my +'SUM')
            calc = rain + melt
            calc.save(path2+newdn)
            print(newdn)
        
calc_avail_water(path, path2)

"""
Explanation: Calculations
Available Water
Calculates monthly available water:
    $$Available\;water = Rain + Snowmelt$$
<br>
End of explanation
"""


from arcpy import Raster

monthRange = [1,12] 
yearRange = [2004,2014]

path = "U:/GWP/Groundwater/Projects/BCM/Data/"
field_cap = Raster(path + "Soil.gdb/fieldCap")
wilt_point = Raster(path + "Soil.gdb/WiltPoint")
T_soil_water = Raster(path + "Soil.gdb/Tsoilwater")
geol_k = Raster(path + "Soil.gdb/Geol_K")
#geol_k = Raster(path + "Soil.gdb/BMC_K")
avail_water = "U:/GWP/Groundwater/Projects/BCM/Data/AvailableWater2.gdb/"
pet = path+'MODIS16.gdb/'

results = "U:/GWP/Groundwater/Projects/BCM/Data/Results.gdb/"
    
def UBM_calc(results, field_cap, wilt_point, T_soil_water, geol_k, avail_water, pet, months = '', years = ''):
    
    arcpy.env.workspace = path
    arcpy.env.overwriteOutput = True

    if months == '':
        months = [1,12] 
    if years == '':
        years = [2004,2014]

    av_soil_water = field_cap
    for y in range(years[0],years[1]+1): #set years converted here
        for m in range(months[0],months[1]+1): #set months converted here
            my = str(y) + str(m).zfill(2)
            av_wtr = Raster(avail_water+'AVWT'+ my)
            pet_rast = Raster(pet + 'PET'+ my) 

            av_soil_water = av_wtr + av_soil_water

            #Eq 1
            av_recharge1 = "T_soil_water - field_cap"
            recharge1 = "Con(eval(av_recharge1) > geol_k, geol_k, eval(av_recharge1))"
            runoff1 = "(av_soil_water - T_soil_water) + Con(eval(av_recharge1) > geol_k, eval(av_recharge1) - geol_k, 0)"
            #Eq2
            av_recharge2 = "av_soil_water - field_cap"
            recharge2 = "Con(eval(av_recharge2) > geol_k, geol_k, eval(av_recharge2))"
            runoff2 = "Con(eval(av_recharge2) > geol_k, eval(av_recharge2) - geol_k, 0)"
            #Eq3 recharge3 = 0 runoff3 = 0 aet = pet_rast
            #Eq4 recharge3 = 0 runoff3 = 0 aet = 0

            #Order of if/then is Eq 1, Eq 4, Eq 2, Eq 3
            recharge = Con(av_soil_water > T_soil_water, eval(recharge1), 
                      Con(av_soil_water < wilt_point, 0,
                         Con((av_soil_water < T_soil_water) & (av_soil_water > field_cap), eval(recharge2),
                            Con((av_soil_water > wilt_point) & (av_soil_water < field_cap), 0, 0))))
            recharge.save(results + 'rec' + my)

            runoff = Con(av_soil_water > T_soil_water, eval(runoff1), 
                      Con(av_soil_water < wilt_point, 0,
                         Con((av_soil_water < T_soil_water) & (av_soil_water > field_cap), eval(runoff2),
                            Con((av_soil_water > wilt_point) & (av_soil_water < field_cap), 0, 0))))
            runoff.save(results + 'run' + my)

            aet = Con(av_soil_water > T_soil_water, pet_rast, 
                      Con(av_soil_water < wilt_point, 0,
                         Con((av_soil_water < T_soil_water) & (av_soil_water > field_cap), pet_rast,
                            Con((av_soil_water > wilt_point) & (av_soil_water < field_cap), pet_rast,0))))
            aet.save(results + 'aet' + my)

            av_soil_water = av_soil_water - runoff - recharge - aet
            av_soil_water = Con(av_soil_water > 0.0, av_soil_water, 0.0)
            av_soil_water.save(results + 'asw' + my)
            print(my)

from arcpy import env  
env.overwriteOutput = True  
env.workspace = "C:/Users/PAULINKENBRANDT/AppData/Roaming/ESRI/Desktop10.4/ArcCatalog/AGRC_SGID.sde"

def monthly_to_yearly(path, code, yearRange='', statistics='SUM'):
    if yearRange=='':
        yearRange = [2004,2014]
    arcpy.env.workspace = path
    arcpy.env.overwriteOutput = True
    for y in range(yearRange[0],yearRange[1]+1): #set years converted here
        ylist = arcpy.ListRasters(code+str(y)+"*") #pick all files from raw data folder of a data type
        calc = CellStatistics(ylist, statistics_type = statistics, ignore_nodata="DATA")
        outnm = 'y'+code+str(y)
        desc = arcpy.Describe(calc)
        print(outnm)
        calc.save(outnm)


def summarize(path, code, statistics='MEAN'):
    arcpy.env.workspace = path
    arcpy.env.overwriteOutput = True

    rlist = arcpy.ListRasters(code+"*") #pick all files from raw data folder of a data type
    # arcpy sa functions that summarize the daily data to monthly data
    calc = CellStatistics(rlist, statistics_type = statistics, ignore_nodata="DATA")
    outnm = code+"_"+statistics
    calc.save(outnm)
    print(outnm)
    

path = 'H:/GIS/Results.gdb'
code = 'rec'
monthly_to_yearly(path,code)
#summarize(path,code)

code = 'yrec'
summarize(path,code)

path = 'H:/GIS/Results.gdb'
code = 'run'
monthly_to_yearly(path,code)
code = 'yrun'
summarize(path,code)

"""
Explanation: Available Soil Water
<p>The calculation of excess water provides the water that is available for watershed hydrology. Available water occupies the soil profile, where water will become actual evapotranspiration,and may result in runoff or recharge, depending on the permeability of the underlying bedrock. Total soil‐water storage is calculated as porosity multiplied by soil depth. Field capacity (soil water volume at ‐0.3 megapascals [MPa]) is the soil water volume below which drainage is negligible, and wilting point (soil water volume at ‐1.5 MPa) is the soil water volume below which actual evapotranspiration does not occur (Hillel 1980).</p>
<p>Once available water is calculated, it may exceed total soil storage and become runoff, or it may be less than total soil storage but greater than field capacity and become recharge. Anything less than field capacity will be calculated as actual evapotranspiration at the rate of PET for that month until it reaches wilting point. When soil water is less than total soil storage and greater than field capacity, soil water greater than field capacity equals recharge. If recharge is greater than bedrock permeability (K), then recharge = K and excess becomes runoff, else it will recharge at K until field capacity. Runoff and recharge combine to calculate basin discharge, and actual evapotranspiration is subtracted from PET to calculate climate water deficit.</p>

$$Available\;Soil\;Water_{t} = Available\;Water_{t} + Available\;Soil\;Water_{t-1}$$
$$Available\;Soil\;Water_{t} = Available\;Soil\;Water_{t-1} - Runoff_{t} - Recharge_{t} - AET_{t}$$
* IF Available Soil Water is greater than Total Soil Water then equation 1<br>
* IF Available Soil Water is between Total Soil Water and FC Soil Water then equation 2<br>
* IF Available Soil Water is between FC Soil Water and WLT Soil Water then equation 3<br>
* IF Available Soil Water is less than WLT Soil Water then equation 4<br>
* Current Month Soil Water becomes new Previous Month Soil Water
<strong>Equation 1</strong>
* Runoff = ((Previous Soil Water + Available Water) - Total Soil Water) + (If recharge > geologic K then (Recharge-geologic K)) 
* Recharge = (Total Soil Water - FC Soil Water) 
* AET = PET
<strong>Equation 2</strong>
* Runoff = (If recharge > geologic K then (Recharge-geologic K)) 
* Recharge = Soil Water - FC Soil Water
* AET = PET
<strong>Equation 3</strong>
* Runoff = 0
* Recharge = 0
* AET = PET
<strong>Equation 4</strong>
* Runoff = 0
* Recharge = 0
* AET = 0
http://gis.stackexchange.com/questions/95027/multiple-conditional-statement-with-in-a-con-tool
End of explanation
"""


