/
Raplee_Kaibab Jobs

Raplee_Kaibab Jobs

Aug 29, 2016

Background

Weathering limited hillslopes exhibit topographic asymmetry in both the size and proportion of cliff abundances along opposing aspects, i.e. south and west facing hillslopes versus north and east facing hillslopes. The mechanisms which lead to these formations are thought to be driven by variations in the net energy balance of the landscape. The energy balance is influenced by differences in hillslope angle, aspect of exposure, and location relative to other topographic features (i.e. daily topographic shading) which influence the timing and quantity of solar insolation and the average air and surface temperatures.

GIS algorithms for calculating solar insolation are resource intensive and time consuming. To alleviate the bottleneck in computation time calculations were performed using a distributed computing workflow which uses the Cooperative Computing Lab’s (CCL) Makeflow (Albrect et al. 2012) and Work Queue (Bui et al. 2011) on the San Diego Super Computing (SDSC) Comet.

Hillslopes in the northern hemisphere tend to have the highest vegetation cover on north and east aspects, and the least on south and west aspects. We hypothesized that this pattern is partly a result of potential evapotranspiration being a function of the product of solar insolation and a nonlinear function of temperature. We evaluate this hypothesis by quantifying potential evapotranspiration in a model framework that resolves the diurnal cycles of both solar insolation and temperature. 

Here we propose a modification to the Hamon equation (Hamon 1961, Haith and Shoenaker 1987), 

PET = 2.1Ht2es / Tt+273.2

where

PET is the potential evaporation/transpiration [mm day-1]

Ht is the average number of hours daylight day-1 during the month of day t

es  saturated vapor pressure at temperature T [kPa], calculated from the es_min  and es_max

Tt is the mean temperature on day t [C] - which we influence using the solar insolation [Si] ratio against mean monthly temperature. Si is the ratio of total_sun_day / flat_sun_day

In the proposed model we take the ratio Si and modify surface temperature based on the overall proportion of daily solar irradiation relative to a flat unshaded surface at the reference elevation (calculated from GRASS r.sun). we also use the total number of hours of sunlight calculated by r.sun in the Ht  term. 

To normalize the total (global) insolation received at the surface and the EarthEnv MODIS derived monthly cloudiness layer (Wilson and Jetz 2016) to correct for the clear-sky radiation model. 

Global radiation * (1 - Cloudiness/2).

Methods

 

GDAL warping and coarsening. Converted the projection to WGS84 UTM Zone 12N, adjusted the pixel size to 16m and 32m.

gdalwarp -overwrite -s_srs "+proj=lcc +lat_1=25 +lat_2=60 +lat_0=42.5 +lon_0=-100 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs" -t_srs EPSG:32612 -r near -tr 16 16 -multi -of GTiff F:\raplee_kaibab\raplee\hand_run\raplee_annual_pet_trad.tif F:/raplee_kaibab/raplee/PET/annual_pet.tif

 

 

DEM generation

The Raplee Ridge lidar were flown by the National Center for Airborne Laser Mapping (NCALM); those data are hosted by the OpenTopography.org data service. A 2 meter (m) DEM was generated by OpenTopo using Martin Isenberg's triangular irregular network (TIN) tin2dem method.

The Kaibab Plateau DEM data were flown by Watershed Sciences and provided by the United States Forest Service Region 3 office. A two meter DEM was generated by interpolation of the original 1m vendor derived bare earth surface model. 

Solar irradiation

To calculate monthly solar irradiation and average monthly temperatures we used an open-source GIS (GRASS v6.4 2015) solar insolation algorithm called r.sun (Hofierka et al. 2007). R.sun calculates daily global irradiation (direct beam, diffuse, and reflected) for clear sky conditions. R.sun also incorporates topographic shading from nearby features. We summed the daily layers, calculated using the default settings for surface albedo (0.2) and Linke atmospheric turbidity (3.0) with a 15 minute sun angle interval to monthly solar irradiation maps. The integration of daily values across an entire month reduces the angular step-spacing between 15 minute intervals as the day lengths change over the course of the month. 

Example of a single day (Julian day 289, October 16) versus monthly average for the month of October. The 15 minute interval for calculating sun hours per day is obvious in the Julian date, but is less clear in the Global radiation surface (although it is still present)

 Daily (Julian Day 289)Monthly (October)
Sun Hours (hours day)
Global Radiation (w m2)

 

For the climate variable we use the 1km DAYMET data product (Thornton et al. 2014) for monthly temperature (minimum and maximum) averaged for the years 1980-2014. The DAYMET data are interpolated using a bicubic spline to match the resolution of the 2m DEM.

Example of two predominantly facing southwest versus southeast aspect cliffs along the eastern edge of the Kaibab Plateau.

 SW Aspect (Top Left) vs SE Aspect (Bottom Left)
Max Temperature (C)
Effective Precipitation (cm month)
  
PET (mm day-1)

 

We processed the Raplee DEM (38 Mb) on Comet using the effective energy and mass transfer (EEMT) workflow available on GitHub: https://github.com/tyson-swetnam/eemt.

Scripts

After the jobs completed on Comet we wrote two shell scripts to calculate the high resolution temperature models in QGIS.

Import and upscale the rasters:

import_eemt_input_raplee.sh

import_eemt_input_kaibab.sh

Calculate input parameters for EEMT, including topographically modified temperature:

calculate_local_temps.sh

Calculating local temperature using solar irradiation:

rsun_modified_hamon_equation.sh

QGIS Python Code for calculating modified Hamon's Equation:

#######################################
###This script imports the necessary layers for calculating Hamon Equation PET
from osgeo import gdal
from osgeo import osr
from PyQt4.QtCore import QFileInfo, QSettings
from qgis.core import QgsRasterLayer, QgsCoordinateReferenceSystem
from qgis.analysis import QgsRasterCalculator, QgsRasterCalculatorEntry
import processing,math,os,sys

## User must set working directory and the filename of the DEM here before executing run
###Set the working directory for the input data
input_Dir = 'F:/raplee_kaibab/kaibab/'
###Set the output directory for the products
output_Dir = 'F:/raplee_kaibab/kaibab/hand_run/'
os.mkdir(output_Dir, 0755)
dem_filename = 'kaibab_converted.tif'

#########################################
##GDAL import DEM to read the projection
###dem=raster
###e.g. dem = 'dem_10m.tif'
dem = gdal.Open(os.path.join(input_Dir,dem_filename))
gt = dem.GetGeoTransform()
xorigin, xcellsize, xrotation, yorigion, yrotation, ycellsize = gt

##Get and Print Projection Information
projInfo=dem .GetProjection()
print projInfo

###Establish number of days per month (for daily calculations to montly)
month_days = {'jan': 31,'feb': 28, 'mar': 31, 'apr': 30, 'may': 31, 'jun': 30, 'jul':31, 'aug': 31, 'sep': 30, 'oct': 31, 'nov': 30, 'dec': 31}
demLayer  = QgsRasterLayer(os.path.join(input_Dir,dem_filename), "dem")
if not demLayer.isValid(): print "DEM Layer failed to load!"

##Load raster layer into canvas
QgsMapLayerRegistry.instance().addMapLayer(demLayer)
 
###Print Extent of the DEM layer (for use with raster calculator)
print demLayer.extent().xMinimum(), demLayer.extent().xMaximum()
print demLayer.extent().yMinimum(), demLayer.extent().yMaximum()
print demLayer.width(), demLayer.height()
name = demLayer.name()
print name
bands = demLayer.bandCount()
print bands

########Local Topographically Corrected Temperature (modified by using the global solar radiation)
print "Calculating locally corrected min and max temperature (maximum temp modified by solar angle, min temp assumed at night)"
jan_S_i = processing.runalg("saga:rastercalculator",os.path.join(input_Dir,'global/monthly/flat_total_sun_jan_sum.tif'),os.path.join(input_Dir,'global/monthly/total_sun_jan_sum.tif'),"b/a",False,8,os.path.join(output_Dir,'jan_S_i.sdat'))
src_ds = gdal.Open(os.path.join(output_Dir,'jan_S_i.sdat'))
##Open output format driver, see gdal_translate --formats for list
format = "GTiff"
driver = gdal.GetDriverByName( format )

##Output to new format
dst_ds = driver.CreateCopy(os.path.join(output_Dir,'jan_S_i.tif'), src_ds, 0 )

##Properly close the datasets to flush to disk
dst_ds = None
src_ds = None
jan_S_i = os.path.join(output_Dir,'jan_S_i.tif')
jan_S_iLayer  = QgsRasterLayer(jan_S_i, "jan_S_i")
if not jan_S_iLayer.isValid(): print "jan_S_i Layer failed to load!"
####Load raster layer into canvas
QgsMapLayerRegistry.instance().addMapLayer(jan_S_iLayer)
######Topographically modified minimum temperatures are at night - values do not include solar radiation e.g. tmin_local = tmin_topo
######Topographically modified max temperature
jan_tmax_topo = processing.runalg("saga:rastercalculator",os.path.join(output_Dir,'jan_tmax_local.tif'),os.path.join(output_Dir,'jan_S_i.tif'),"a+(b-(1/b))",False,8,os.path.join(output_Dir,'jan_tmax_topo.sdat'))
src_ds = gdal.Open(os.path.join(output_Dir,'jan_tmax_topo.sdat'))
##Open output format driver, see gdal_translate --formats for list
format = "GTiff"
driver = gdal.GetDriverByName( format )

##Output to new format
dst_ds = driver.CreateCopy(os.path.join(output_Dir,'jan_tmax_topo.tif'), src_ds, 0 )
##Properly close the datasets to flush to disk
dst_ds = None
src_ds = None
jan_tmax_topo = os.path.join(output_Dir,'jan_tmax_topo.tif')
jan_tmax_topoLayer  = QgsRasterLayer(jan_tmax_topo, "jan_tmax_topo")
if not jan_tmax_topoLayer.isValid(): print "jan_tmax_topo Layer failed to load!"

####Load raster layer into canvas
QgsMapLayerRegistry.instance().addMapLayer(jan_tmax_topoLayer)

######Topographically modified mean temperature
jan_tmean_topo = processing.runalg("saga:rastercalculator",os.path.join(output_Dir,'jan_tmax_topo.tif'),os.path.join(output_Dir,'jan_tmin_local.tif'),"(a+b)/2",False,8,os.path.join(output_Dir,'jan_tmean_topo.sdat'))
src_ds = gdal.Open(os.path.join(output_Dir,'jan_tmean_topo.sdat'))
##Open output format driver, see gdal_translate --formats for list
format = "GTiff"
driver = gdal.GetDriverByName( format )
##Output to new format
dst_ds = driver.CreateCopy(os.path.join(output_Dir,'jan_tmean_topo.tif'), src_ds, 0 )
##Properly close the datasets to flush to disk
dst_ds = None
src_ds = None
jan_tmean_topo = os.path.join(output_Dir,'jan_tmean_topo.tif')
jan_tmean_topoLayer  = QgsRasterLayer(jan_tmean_topo, "jan_tmean_topo")
if not jan_tmean_topoLayer.isValid(): print "jan_tmean_topo Layer failed to load!"
####Load raster layer into canvas
QgsMapLayerRegistry.instance().addMapLayer(jan_tmean_topoLayer)

########Begin calculating the input parameters for PET
print "Calculating local Potential Evaporation-Transpiration (PET) for EEMT-Trad using Hamon's Equation"
######Calculate es_tmax
jan_es_tmax = processing.runalg("saga:rastercalculator",os.path.join(output_Dir,'jan_tmax_topo.tif'),None,"0.6108*exp((17.27*a)/(a+237.3))",False,8,os.path.join(output_Dir,'jan_es_tmax.sdat'))
src_ds = gdal.Open(os.path.join(output_Dir,'jan_es_tmax.sdat'))
##Open output format driver, see gdal_translate --formats for list
format = "GTiff"
driver = gdal.GetDriverByName( format )
##Output to new format
dst_ds = driver.CreateCopy(os.path.join(output_Dir,'jan_es_tmax.tif'), src_ds, 0 )
##Properly close the datasets to flush to disk
dst_ds = None
src_ds = None
jan_estmax = os.path.join(output_Dir,'jan_es_tmax.tif')
jan_estmaxLayer  = QgsRasterLayer(jan_estmax, "jan_estmax")
if not jan_estmaxLayer.isValid(): print "jan_es_tmax Layer failed to load!"
####Load raster layer into canvas
QgsMapLayerRegistry.instance().addMapLayer(jan_estmaxLayer)
######Calculate es_tmin
jan_es_tmin = processing.runalg("saga:rastercalculator",os.path.join(output_Dir,'jan_tmin_local.tif'),None,"0.6108*exp((17.27*a)/(a+237.3))",False,8,os.path.join(output_Dir,'jan_es_tmin.sdat'))
src_ds = gdal.Open(os.path.join(output_Dir,'jan_es_tmin.sdat'))
##Open output format driver, see gdal_translate --formats for list
format = "GTiff"
driver = gdal.GetDriverByName( format )
##Output to new format
dst_ds = driver.CreateCopy(os.path.join(output_Dir,'jan_es_tmin.tif'), src_ds, 0 )
##Properly close the datasets to flush to disk
dst_ds = None
src_ds = None
jan_estmin = os.path.join(output_Dir,'jan_es_tmin.tif')
jan_estminLayer  = QgsRasterLayer(jan_estmin, "jan_estmin")
if not jan_estminLayer.isValid(): print "jan_es_tmin Layer failed to load!"

####Load raster layer into canvas
QgsMapLayerRegistry.instance().addMapLayer(jan_estminLayer)
######Calculate the mean saturated vapor pressure (kPA)
####Mask pixels with lower than zero mean monthly temperatures for e_s
## When there are more than 2 inputs you need to generate a separate input_params input
input_params =[os.path.join(output_Dir,'jan_es_tmin.sdat'),os.path.join(output_Dir,'jan_es_tmax.sdat')]
jan_e_s = processing.runalg("saga:rastercalculator",os.path.join(output_Dir,'jan_tmean_topo.sdat'),input_params,"ifelse(lt(a,0),0,(b+c)/2)",False,8,os.path.join(output_Dir,'jan_e_s.sdat'))
src_ds = gdal.Open(os.path.join(output_Dir,'jan_e_s.sdat'))
##Open output format driver, see gdal_translate --formats for list
format = "GTiff"
driver = gdal.GetDriverByName( format )
##Output to new format
dst_ds = driver.CreateCopy(os.path.join(output_Dir,'jan_e_s.tif'), src_ds, 0 )
##Properly close the datasets to flush to disk
dst_ds = None
src_ds = None
jan_e_s = os.path.join(output_Dir,'jan_e_s.tif')
jan_esLayer  = QgsRasterLayer(jan_e_s, "jan_e_s")
if not jan_esLayer.isValid(): print "jan_e_s Layer failed to load!"
####Load raster layer into canvas
QgsMapLayerRegistry.instance().addMapLayer(jan_esLayer)

########PET (Traditional) calculation WARNING: Units are in mm^-1 day^-1 
########see https://dl.sciencesocieties.org/publications/vzj/supplements/14/vzj2014.07.0102-supplement.pdf
input_params=[os.path.join(output_Dir,'jan_e_s.tif'),os.path.join(output_Dir,'jan_tmean_topo.tif')]
jan_PET_trad = processing.runalg("saga:rastercalculator",os.path.join(input_Dir,'insol/monthly/hours_sun_jan_average.tif'),input_params,"ifelse(eq(b,0),0,(2.1*(a^2)*b)/(c+273.2))",False,8,os.path.join(output_Dir,'jan_PET_trad.sdat'))
src_ds = gdal.Open(os.path.join(output_Dir,'jan_PET_trad.sdat'))
##Open output format driver, see gdal_translate --formats for list
format = "GTiff"
driver = gdal.GetDriverByName( format )
##Output to new format
dst_ds = driver.CreateCopy(os.path.join(output_Dir,'jan_PET_trad.tif'), src_ds, 0 )
##Properly close the datasets to flush to disk
dst_ds = None
src_ds = None
jan_PET_trad = os.path.join(output_Dir,'jan_PET_trad.tif')
jan_PETtradLayer  = QgsRasterLayer(jan_PET_trad, "jan_PET_trad")
if not jan_PETtradLayer.isValid(): print "jan_PET_trad Layer failed to load!"
####Load raster layer into canvas
QgsMapLayerRegistry.instance().addMapLayer(jan_PETtradLayer)

Results

Comet Jobs

The job ran smoothly, taking ~ 5.27 hours in wall time; we was charged 827 SU.

The solar insolation output of the Raplee EEMT is ~44.5 Gb

Job time results:

real    316m9.945s
user    0m25.089s
sys     1m36.359s

Next, we ran the Kaibab DEM (98 Mb); there were several errors running the job on the Comet login. These are apparently caused by Comet's limit on the number of jobs it can have open at one time. After talking with Mats we ran the jobs on the xd-login (head node) using the Comet resources.

The job ran smoothly and completed. Job time results:

real    665m45.825s
user    0m46.228s
sys     3m15.032s

The solar insolation output of the Kaibab EEMT is 71 Mb. The output of Kaibab is ~96 Gb. 

Acknowledgements 

This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1053575. We thank Mats Ryge and Yan Liu for their assistance with the EEMT workflow, which was made possible through the XSEDE Extended Collaborative Support Service (ECSS) program.

The Kaibab lidar were funded by the Forest Service Southwestern Region to support landscape restoration and northern Goshawk research on the North Kaibab Plateau.  Lidar QA/QC was by the Southwestern Regional Office Remote Sensing Unit.

The Raplee Ridge lidar data acquisition and processing were completed by the National Center for Airborne Laser Mapping (NCALM - http://www.ncalm.org). NCALM funding provided by NSF's Division of Earth Sciences, Instrumentation and Facilities Program. EAR-1043051.

References

Albrecht, M., P. Donnelly, P. Bui, and D. Thain. 2012.Makeflow: A Portable Abstraction for Data IntensiveComputing on Clusters, Clouds, and Grids. In Workshop onScalable Workflow Enactment Engines and Technologies(SWEET) at ACM SIGMOD.

Bui P., D. Rajan, B. Abdul-Wahid, J. Izaguirre, and D.Thain. 2011. Work Queue + Python: A Framework ForScalable Scientific Ensemble Applications. In Workshop onPython for High Performance and Scientific Computing(PyHPC) at the ACM/IEEE International Conference forHigh Performance Computing, Networking, Storage, andAnalysis (Supercomputing).

GRASS Development Team, 2015. Geographic ResourcesAnalysis Support System (GRASS) Software, Version 6.4. Open Source Geospatial Foundation.http://grass.osgeo.org

Haith, D.A., and L.L. Shoenaker. 1987. Generalized watershed loading functions for stream-flow nutrients. Water Resour. Bull. 23:471–478. doi:10.1111/j.1752-1688.1987.tb00825.x

Hamon W.R. 1961. Estimating potential evapotranspiration. Proc. Am. Soc. Civ. Eng. 87:107–120.

Hofierka, J., T. Huld, T. Cebecauer, and M. Šúri, 2007. Opensource solar radiation tools for environmental and renewableenergy applications. Environmental Software Systems, 448.(doi: 10.13140/2.1.3773.6001

Thornton, P.E., M.M. Thornton, B.W. Mayer, N. Wilhelmi,Y.Wei, R. Devarakonda, and R.B. Cook. 2014. Daymet:Daily Surface Weather Data on a 1-km Grid for NorthAmerica, Version 2. Data set. Available on-line[http://daac.ornl.gov] from Oak Ridge National LaboratoryDistributed Active Archive Center, Oak Ridge, Tennessee,USA.

John Towns, Timothy Cockerill, Maytal Dahan, Ian Foster, Kelly Gaither, Andrew Grimshaw, Victor Hazlewood, Scott Lathrop, Dave Lifka, Gregory D. Peterson, Ralph Roskies, J. Ray Scott, Nancy Wilkins-Diehr, "XSEDE: Accelerating Scientific Discovery", Computing in Science & Engineering, vol.16, no. 5, pp. 62-74, Sept.-Oct. 2014, doi:10.1109/MCSE.2014.80

Wilkins-Diehr, N and S Sanielevici, J Alameda, J Cazes, L Crosby, M Pierce, R Roskies. High Performance Computer Applications 6th International Conference, ISUM 2015, Mexico City, Mexico, March 9-13, 2015, Revised Selected Papers Gitler, Isidoro, Klapp, Jaime (Eds.) Springer International Publishing. ISBN 978-3-319-32243-8, 3-13, 2016. 10.1007/978-3-319-32243-8.

Wilson AM, Jetz W (2016) Remotely Sensed High-Resolution Global Cloud Dynamics for Predicting Ecosystem and Biodiversity Distributions. PLoS Biol 14(3): e1002415. doi:10.1371/journal. pbio.1002415