Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Updated python script #3

Open
wants to merge 2 commits into
base: master
Choose a base branch
from
Open

Updated python script #3

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
6b0098a
Modified script so it could be used as importable module
nina-aragon Aug 21, 2020
a871871
rm extra files
nina-aragon Aug 21, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
106 changes: 49 additions & 57 deletions Depth2SWE.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,9 @@
"metadata": {},
"source": [
"# SWE Calculator \n",
"\n",
"This is a python implementation of a matlab script by D. Hill. See full details at: Hill, D.F., Burakowski, E.A., Crumley, R.L., Keon, J., Hu, J.M., Arendt, A.A., Wikstrom Jones, K., Wolken G.J., 2019, \"Converting snow depth to snow water equivalent using climatological variables,\" The Cryosphere, v.13, pp.1767-1784 https://doi.org/10.5194/tc-13-1767-2019\n",
"\n",
"This script implements a power law regression in order to produce a value of snow water equivalent (SWE) based on various parameters associated with a snow depth measurement. The input variables can be scalars (single measurement) or vectors (batch measurements). The following information on units, etc., is critical. To use this script, edit inputs with your data.\n",
"\n",
"Christina Aragon, Oregon State University \n",
Expand Down Expand Up @@ -43,64 +46,53 @@
"from scipy.interpolate import interp2d\n",
"from datetime import date\n",
"\n",
"#***Edit input data***. Input variables can be a list or array. \n",
"#Below is an example of manually entered inputs, or data can be imported from .txt, .csv, etc. \n",
"Y = [2018,2018,2018]\n",
"M = [1,11,3]\n",
"D = [1,1,1]\n",
"H = [30,40,50]\n",
"LAT = [43.5,43.5,43.5]\n",
"LON = [-110.8,-110.8,-110.8]\n",
"\n",
"#Build lat/lon array \n",
"ncols = 7300\n",
"nrows = 2839\n",
"xll = -168.00051894775\n",
"yll = 30.002598288104\n",
"clsz = 0.014795907586\n",
"#Longitudes\n",
"ln = np.arange(xll, xll+ncols*clsz, clsz)\n",
"#Latitudes\n",
"lt = np.arange(yll, yll+nrows*clsz, clsz)\n",
"la = np.flipud(lt)\n",
"\n",
"#Import temperature difference data\n",
"geo = gdal.Open('td_final.txt')\n",
"td = geo.ReadAsArray()\n",
"#Import winter precipitation data\n",
"geo = gdal.Open('ppt_wt_final.txt')\n",
"pptwt = geo.ReadAsArray()\n",
" \n",
"def swe_calc(Y,M,D,H,LAT,LON):\n",
" #Interpolate temp to input\n",
" f_td = interp2d(ln, la, td, kind='cubic')\n",
" TD = f_td(LON, LAT)\n",
" #Interpolate temp to input\n",
" f_ppt = interp2d(ln, la, pptwt, kind='cubic')\n",
" PPTWT = f_ppt(LON, LAT)\n",
" #Determine day of year\n",
" DOY = date.toordinal(date(Y,M,D))-date.toordinal(date(Y,9,30))\n",
" if DOY < 0:\n",
" DOY = DOY+365\n",
" #Apply regression equation \n",
" a = [0.0533,0.948,0.1701,-0.1314,0.2922] #accumulation phase\n",
" b = [0.0481,1.0395,0.1699,-0.0461,0.1804]; #ablation phase\n",
" SWE = a[0]*H**a[1]*PPTWT**a[2]*TD**a[3]*DOY**a[4]* \\\n",
" (-np.tanh(.01*(DOY-180))+1)/2 + b[0]*H**b[1]* \\\n",
" PPTWT**b[2]*TD**b[3]*DOY**b[4] * (np.tanh(.01*(DOY-180))+1)/2;\n",
" return SWE,DOY\n",
" #Build lat/lon array \n",
" ncols = 7300\n",
" nrows = 2839\n",
" xll = -168.00051894775\n",
" yll = 30.002598288104\n",
" clsz = 0.014795907586\n",
" #Longitudes\n",
" ln = np.arange(xll, xll+ncols*clsz, clsz)\n",
" #Latitudes\n",
" lt = np.arange(yll, yll+nrows*clsz, clsz)\n",
" la = np.flipud(lt)\n",
"\n",
" #Import temperature difference data\n",
" geo = gdal.Open('td_final.txt')\n",
" td = geo.ReadAsArray()\n",
" #Import winter precipitation data\n",
" geo = gdal.Open('ppt_wt_final.txt')\n",
" pptwt = geo.ReadAsArray()\n",
"\n",
"#Create output arrays \n",
"SWE = np.zeros(len(H))*np.nan\n",
"DOY = np.zeros(len(H))*np.nan\n",
"for i in np.arange(len(H)):\n",
" swe, doy = swe_calc(Y[i],M[i],D[i],H[i],LAT[i],LON[i])\n",
" SWE[i] = swe\n",
" DOY[i] = doy\n",
" def calc(Y,M,D,H,LAT,LON):\n",
" #Interpolate temp to input\n",
" f_td = interp2d(ln, la, td, kind='cubic')\n",
" TD = f_td(LON, LAT)\n",
" #Interpolate temp to input\n",
" f_ppt = interp2d(ln, la, pptwt, kind='cubic')\n",
" PPTWT = f_ppt(LON, LAT)\n",
" #Determine day of year\n",
" DOY = date.toordinal(date(Y,M,D))-date.toordinal(date(Y,9,30))\n",
" if DOY < 0:\n",
" DOY = DOY+365\n",
" #Apply regression equation \n",
" a = [0.0533,0.948,0.1701,-0.1314,0.2922] #accumulation phase\n",
" b = [0.0481,1.0395,0.1699,-0.0461,0.1804]; #ablation phase\n",
" SWE = a[0]*H**a[1]*PPTWT**a[2]*TD**a[3]*DOY**a[4]*(-np.tanh(.01*\\\n",
" (DOY-180))+1)/2 + b[0]*H**b[1]*PPTWT**b[2]*TD**b[3]*DOY**b[4]*\\\n",
" (np.tanh(.01*(DOY-180))+1)/2;\n",
" return SWE,DOY\n",
"\n",
"#Optional: print output SWE and DOY \n",
"#print(SWE)\n",
"#print(DOY)"
" #Create output arrays \n",
" SWE = np.zeros(len(H))*np.nan\n",
" DOY = np.zeros(len(H))*np.nan\n",
" for i in np.arange(len(H)):\n",
" swe, doy = calc(Y[i],M[i],D[i],H[i],LAT[i],LON[i])\n",
" SWE[i] = swe\n",
" DOY[i] = int(doy)\n",
" return SWE, DOY"
]
}
],
Expand All @@ -120,9 +112,9 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.0"
"version": "3.8.3"
}
},
"nbformat": 4,
"nbformat_minor": 2
"nbformat_minor": 4
}
118 changes: 56 additions & 62 deletions Depth2SWE.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,19 +1,21 @@
#!/usr/bin/env python
# coding: utf-8

#################### SWE Calculator ####################
# # SWE Calculator
#
# This is a python implementation of a matlab script by D. Hill. See full details at: Hill, D.F., Burakowski, E.A., Crumley, R.L., Keon, J., Hu, J.M., Arendt, A.A., Wikstrom Jones, K., Wolken G.J., 2019, "Converting snow depth to snow water equivalent using climatological variables," The Cryosphere, v.13, pp.1767-1784 https://doi.org/10.5194/tc-13-1767-2019
#
# This script implements a power law regression in order to produce a value of snow water equivalent (SWE) based on various parameters associated with a snow depth measurement. The input variables can be scalars (single measurement) or vectors (batch measurements). The following information on units, etc., is critical. To use this script, edit inputs with your data.

#
# Christina Aragon, Oregon State University
# December 2019
#######################################################################

#################### Data ####################

# Gridded winter precipitation (PPTWT) and temperature difference (TD) data is used in this script. See https://adaptwest.databasin.org/pages/adaptwest-climatena for info on data. See http://www.cacpd.org.s3-website-us-west-2.amazonaws.com/climate_normals/NA_ReadMe.txt for specific metadata.
#
# ## Data
# Gridded winter precipitation (PPTWT) and temperature difference (TD) data is used in this script. See https://adaptwest.databasin.org/pages/adaptwest-climatena for info on data. See http://www.cacpd.org.s3-website-us-west-2.amazonaws.com/climate_normals/NA_ReadMe.txt for specific metadata.
#
# The grids that are included with this calculator have been put into geographic coords and they have been clipped to (72>lat>30) and (-168<lon<-60).
#
##################### User Defined Inputs ####################
# ## User Defined Inputs:
#
# #### H: snow depth [mm]
# #### Y: year
Expand All @@ -24,70 +26,62 @@
# #### LON: longitude [degrees]
# Signed - ex: -120 for N America
#
#################### Outputs ####################
# ## Outputs:
#
# #### SWE: snow water equivalent [mm]
# #### DOY: day of water year [Oct 1 = 1]

########################################################################
# In[ ]:


import numpy as np
from osgeo import gdal
from scipy.interpolate import interp2d
from datetime import date

#***Edit input data***. Input variables can be a list or array.
#Below is an example of manually entered inputs, or data can be imported from .txt, .csv, etc.
Y = [2018,2018,2018]
M = [1,11,3]
D = [1,1,1]
H = [30,40,50]
LAT = [43.5,43.5,43.5]
LON = [-110.8,-110.8,-110.8]
def swe_calc(Y,M,D,H,LAT,LON):
#Build lat/lon array
ncols = 7300
nrows = 2839
xll = -168.00051894775
yll = 30.002598288104
clsz = 0.014795907586
#Longitudes
ln = np.arange(xll, xll+ncols*clsz, clsz)
#Latitudes
lt = np.arange(yll, yll+nrows*clsz, clsz)
la = np.flipud(lt)

#Build lat/lon array
ncols = 7300
nrows = 2839
xll = -168.00051894775
yll = 30.002598288104
clsz = 0.014795907586
#Longitudes
ln = np.arange(xll, xll+ncols*clsz, clsz)
#Latitudes
lt = np.arange(yll, yll+nrows*clsz, clsz)
la = np.flipud(lt)
#Import temperature difference data
geo = gdal.Open('td_final.txt')
td = geo.ReadAsArray()
#Import winter precipitation data
geo = gdal.Open('ppt_wt_final.txt')
pptwt = geo.ReadAsArray()

#Import temperature difference data
geo = gdal.Open('td_final.txt')
td = geo.ReadAsArray()
#Import winter precipitation data
geo = gdal.Open('ppt_wt_final.txt')
pptwt = geo.ReadAsArray()

def swe_calc(Y,M,D,H,LAT,LON):
#Interpolate temp to input
f_td = interp2d(ln, la, td, kind='cubic')
TD = f_td(LON, LAT)
#Interpolate temp to input
f_ppt = interp2d(ln, la, pptwt, kind='cubic')
PPTWT = f_ppt(LON, LAT)
#Determine day of year
DOY = date.toordinal(date(Y,M,D))-date.toordinal(date(Y,9,30))
if DOY < 0:
DOY = DOY+365
#Apply regression equation
a = [0.0533,0.948,0.1701,-0.1314,0.2922] #accumulation phase
b = [0.0481,1.0395,0.1699,-0.0461,0.1804]; #ablation phase
SWE = a[0]*H**a[1]*PPTWT**a[2]*TD**a[3]*DOY**a[4]* (-np.tanh(.01*(DOY-180))+1)/2 + b[0]*H**b[1]* PPTWT**b[2]*TD**b[3]*DOY**b[4] * (np.tanh(.01*(DOY-180))+1)/2;
return SWE,DOY
def calc(Y,M,D,H,LAT,LON):
#Interpolate temp to input
f_td = interp2d(ln, la, td, kind='cubic')
TD = f_td(LON, LAT)
#Interpolate temp to input
f_ppt = interp2d(ln, la, pptwt, kind='cubic')
PPTWT = f_ppt(LON, LAT)
#Determine day of year
DOY = date.toordinal(date(Y,M,D))-date.toordinal(date(Y,9,30))
if DOY < 0:
DOY = DOY+365
#Apply regression equation
a = [0.0533,0.948,0.1701,-0.1314,0.2922] #accumulation phase
b = [0.0481,1.0395,0.1699,-0.0461,0.1804]; #ablation phase
SWE = a[0]*H**a[1]*PPTWT**a[2]*TD**a[3]*DOY**a[4]*(-np.tanh(.01* (DOY-180))+1)/2 + b[0]*H**b[1]*PPTWT**b[2]*TD**b[3]*DOY**b[4]* (np.tanh(.01*(DOY-180))+1)/2;
return SWE,DOY

#Create output arrays
SWE = np.zeros(len(H))*np.nan
DOY = np.zeros(len(H))*np.nan
for i in np.arange(len(H)):
swe, doy = swe_calc(Y[i],M[i],D[i],H[i],LAT[i],LON[i])
SWE[i] = swe
DOY[i] = doy
#Create output arrays
SWE = np.zeros(len(H))*np.nan
DOY = np.zeros(len(H))*np.nan
for i in np.arange(len(H)):
swe, doy = calc(Y[i],M[i],D[i],H[i],LAT[i],LON[i])
SWE[i] = swe
DOY[i] = int(doy)
return SWE, DOY

#Optional: print output SWE and DOY
#print(SWE)
#print(DOY)