fix up attribute function

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: mda.lakes
 Type: Package
 Title: Tools for combining models, data, and processing for lakes
-Version: 3.0.2
+Version: 3.0.4
 Date: 2015-12-03
 Author: Luke Winslow, Jordan Read
 Maintainer: Luke Winslow <[email protected]>

NAMESPACE

@@ -1,4 +1,4 @@
-# Generated by roxygen2 (4.1.1): do not edit by hand
+# Generated by roxygen2: do not edit by hand
 
 export(LTER_NLDAS_gapper)
 export(area_light_temp_threshold)
@@ -51,6 +51,7 @@ export(prep_run_chained_glm_kd)
 export(prep_run_glm_kd)
 export(sens_seasonal_site)
 export(set_driver_url)
+export(write_error_log)
 import(GLMr)
 import(glmtools)
 import(lakeattributes)

R/lake_attribute_functions.R

@@ -122,7 +122,7 @@ getResidenceTime	<-	local(
 #' 
 #' 
 #' @export
-getCanopy = function(site_id, default.if.null=FALSE, method="aster"){
+getCanopy = function(site_id, default.if.null=FALSE, method="landcover"){
 
 	if (tolower(method) == 'aster'){
 
@@ -137,28 +137,13 @@ getCanopy = function(site_id, default.if.null=FALSE, method="aster"){
 		return(vals$MEAN_HT)
 
 	}else if (tolower(method) == "landcover"){
-
-		fname = system.file('supporting_files/buffers_land_cover.csv', 
-												package=packageName())
-		d = read.table(fname, header=TRUE, sep=',')
-		#100 urban
-		#110 ag
-		#150 grassland
-		#160 forest 
-		#200 open water
-		#210 wetland
-		lc_types = data.frame(majority_cover=c(200, 160,  100, 150,  110, 210),
-													height=        c(0.5, 11.5, 0.5, 0.65, 0.8, 0.5))
-
-		vals = merge(data.frame(site_id, order=1:length(site_id)),
-								 d, by='site_id', all.x=TRUE)
-
-		heights = merge(vals, lc_types, all.x=TRUE, by='majority_cover')
-
-		#merge does not preserve order, re-order before returning
-		heights = heights[order(heights$order),]
-
-		return(heights$height)
+		data(canopy)
+	  id = site_id
+		d = subset(canopy, `site_id` == id & source == 'nlcd')
+		if(nrow(d) < 1){
+		  return(NA)
+		}
+		return(d$canopy_m)
 
 	}else{
 		stop('Unidentified method ', method, ' for getCanopy [aster, landcover]')
@@ -494,9 +479,9 @@ getCD	<-	function(site_id=NULL, Wstr=NULL, method='Hondzo'){
 #'
 #'
 #'@export
-getWstr	<-	function(WBIC, method='Markfort', canopy=NULL){
+getWstr	<-	function(site_id, method='Markfort', canopy=NULL){
 
-	lkeArea	<-	get_area(WBIC)
+	lkeArea	<-	get_area(site_id)
 
   if(is.na(lkeArea)){
     return(NA)
@@ -506,7 +491,7 @@ getWstr	<-	function(WBIC, method='Markfort', canopy=NULL){
 		# Markfort et al. 2010
 		minWstr	<-	0.0001
 		if (is.null(canopy)){
-			hc	<-	max(c(getCanopy(WBIC),1))
+			hc	<-	max(c(getCanopy(site_id),1))
 		} else {
 			hc	<-	canopy
 		}
@@ -530,14 +515,14 @@ getWstr	<-	function(WBIC, method='Markfort', canopy=NULL){
 		# Markfort et al. 2010
 		minWstr	<-	0.0001
 		if (is.null(canopy)){
-			hc	<-	max(c(getCanopy(WBIC),1))
+			hc	<-	max(c(getCanopy(site_id),1))
 		} else {
 			hc	<-	canopy
 		}
 
 		Xt	<-	50*hc
 
-		perim	<-	getPerim(WBIC)
+		perim	<-	getPerim(site_id)
 		shelArea	<-	perim*hc*12.5 # 25% of Markfort, as sheltering is single direction...
 		shelter	<-	(lkeArea-shelArea)/lkeArea
 		wind.shelter	<-	max(c(shelter,minWstr))

demo/cluster_cont_calibrate_bias.R

@@ -6,12 +6,15 @@
 #local_url = '192.0.2.1:4040'
 
 local_url = paste0((Sys.info()["nodename"]),':4040')
+driver_url = paste0('http://', (Sys.info()["nodename"]),':4041/')
+
 library(parallel)
 
 #lets try 50 to start
-c1 = makePSOCKcluster(paste0('licon', 1:50), manual=TRUE, port=4045)
+c1 = makePSOCKcluster(paste0('licon', 1:50), manual=TRUE, port=4046)
 
 clusterExport(c1, varlist = 'local_url')
+clusterExport(c1, varlist = 'driver_url')
 
 #clusterCall(c1, function(){install.packages('devtools', repos='http://cran.rstudio.com')})
 clusterCall(c1, function(){install.packages('rLakeAnalyzer', repos='http://cran.rstudio.com')})
@@ -24,8 +27,8 @@ clusterCall(c1, function(){library(devtools)})
 #glmr_install     = clusterCall(c1, function(){install.packages('glmtools', repos=c('http://owi.usgs.gov/R', 'http://cran.rstudio.com'))})
 glmr_install     = clusterCall(c1, function(){install_url(paste0('http://', local_url,'/GLMr_3.1.10.tar.gz'))})
 glmtools_install = clusterCall(c1, function(){install_url(paste0('http://', local_url,'/glmtools_0.13.0.tar.gz'))})
-lakeattr_install = clusterCall(c1, function(){install_url(paste0('http://', local_url,'/lakeattributes_0.7.4.tar.gz'))})
-mdalakes_install = clusterCall(c1, function(){install_url(paste0('http://', local_url,'/mda.lakes_3.0.2.tar.gz'))})
+lakeattr_install = clusterCall(c1, function(){install_url(paste0('http://', local_url,'/lakeattributes_0.8.4.tar.gz'))})
+mdalakes_install = clusterCall(c1, function(){install_url(paste0('http://', local_url,'/mda.lakes_3.0.3.tar.gz'))})
 
 library(lakeattributes)
 library(mda.lakes)
@@ -34,7 +37,7 @@ library(glmtools)
 
 ## setup run function
 
-run_cal = function(site_id, years=1979:2012, driver_function=get_driver_path, nml_args=list()){
+run_cal = function(site_id, years=1979:2012, driver_function=get_driver_path, nml_args=list(), datasource=NA){
 
 
   library(lakeattributes)
@@ -51,7 +54,7 @@ run_cal = function(site_id, years=1979:2012, driver_function=get_driver_path, nm
   tryCatch({
 
 
-    run_dir = file.path(fastdir, site_id)
+    run_dir = file.path(fastdir, paste0(site_id, '_', sample.int(1e9, size=1)))
     cat(run_dir, '\n')
     dir.create(run_dir)
 
@@ -71,7 +74,7 @@ run_cal = function(site_id, years=1979:2012, driver_function=get_driver_path, nm
     driver_path = gsub('\\\\', '/', driver_path)
 
 
-    kds = get_kd_best(site_id, years=years)
+    kds = get_kd_best(site_id, years=years, datasource = datasource)
 
     kd_avg = secchi_conv/mean(kds$secchi_avg, na.rm=TRUE)
 
@@ -90,11 +93,14 @@ run_cal = function(site_id, years=1979:2012, driver_function=get_driver_path, nm
     cal.data = resample_to_field(file.path(run_dir, 'output.nc'), file.path(run_dir,'obs.tsv'))
     cal.data$site_id = site_id
 
-    #unlink(run_dir, recursive=TRUE)
+    unlink(run_dir, recursive=TRUE)
 
     return(cal.data)
 
-  }, error=function(e){unlink(run_dir, recursive=TRUE);e})
+  }, error=function(e){
+      unlink(run_dir, recursive=TRUE)
+      return(paste(e$call, e$message, Sys.info()["nodename"]))
+    })
 }
 
 
@@ -103,6 +109,7 @@ driver_fun = function(site_id){
   drivers = driver_nldas_wind_debias(nldas)
   drivers = driver_add_burnin_years(drivers, nyears=2)
   drivers = driver_add_rain(drivers, month=7:9, rain_add=0.5) ##keep the lakes topped off
+  #drivers$time = drivers$time - as.difftime(2, units='days')
   driver_save(drivers)
 }
 
@@ -111,8 +118,14 @@ driver_fun = function(site_id){
 #c1 = c1[unlist(ramdisk)]
 
 to_run = unique(intersect(intersect(wtemp$site_id, zmax$site_id), secchi$site_id))
+#to_run = to_run[to_run %in% read.table('~/managed_lakes_wilma_nhd.tsv', sep='\t', header=TRUE, as.is=TRUE)$id]
+
+run_name = '2016-03-31_cal_NLDAS'
+run_comment = "Full calibration run with all available wtemp, secchi and depth"
+
+clusterCall(c1, function(){library(mda.lakes);set_driver_url(driver_url)})
 
-out = clusterApplyLB(c1, to_run, run_cal, driver_function = driver_fun, nml_args=list(), years=1980:1999)
+out = clusterApplyLB(c1, to_run, run_cal, driver_function = driver_fun, nml_args=list(), years=1980:1999, datasource=NA)
 
 #out = clusterApplyLB(c1, to_run[1:50], run_cal)
 
@@ -121,10 +134,13 @@ sprintf('%i lakes ran\n', sum(unlist(lapply(out, inherits, what='data.frame'))))
 
 ##results
 out_df = do.call('rbind', out[unlist(lapply(out, inherits, what='data.frame'))])
-write.table(out_df, 'c:/WiLMA/results/2016-03-22_NLDAS.tsv', sep='\t', row.names=FALSE)
+write.table(out_df, paste0('c:/WiLMA/results/', run_name, '.tsv'), sep='\t', row.names=FALSE)
 
-rmarkdown::render('demo/document_calibration_overview_template.Rmd', output_file='2016-03-22_NLDAS.pdf', 
-                  output_dir='c:/WiLMA/results/', params=list(out_df=out_df))
+out_err = out[!unlist(lapply(out, inherits, what='data.frame'))]
+save(out_err, file = paste0('c:/WiLMA/results/', run_name, '.Rdata'))
+
+rmarkdown::render('demo/document_calibration_overview_template.Rmd', output_file=paste0(run_name,'.pdf'), 
+                  output_dir='c:/WiLMA/results/', params=list(out_df=out_df, run_message=run_comment))
 
 
 sqrt(mean((out_df$Observed_temp - out_df$Modeled_temp)^2, na.rm=TRUE))
@@ -133,31 +149,39 @@ mean(out_df$Observed_temp - out_df$Modeled_temp, na.rm=TRUE)
 ################################################################################
 ## ACCESS
 ################################################################################
-driver_fun = function(site_id, gcm){
-  drivers = read.csv(get_driver_path(site_id, driver_name = gcm, timestep = 'daily'), header=TRUE)
-  nldas   = read.csv(get_driver_path(site_id, driver_name = 'NLDAS', timestep = 'daily'), header=TRUE)
-  drivers = driver_nldas_debias_airt_sw(drivers, nldas)
+driver_fun = function(site_id){
+  drivers = read.csv(get_driver_path(site_id, driver_name = 'ACCESS', timestep = 'daily'), header=TRUE)
+  #nldas   = read.csv(get_driver_path(site_id, driver_name = 'NLDAS', timestep = 'daily'), header=TRUE)
+  #drivers = driver_nldas_debias_airt_sw(drivers, nldas)
   drivers = driver_add_burnin_years(drivers, nyears=2)
   drivers = driver_add_rain(drivers, month=7:9, rain_add=0.5) ##keep the lakes topped off
   driver_save(drivers)
 }
 
+run_name = '2016-03-28_ACCESS_WILMALakes'
+run_comment = "Look at bias and RMSE for Notaro ACCESS drivers on WiLMA only lakes. Years are only 1980:1999"
+
+clusterCall(c1, function(){library(mda.lakes);set_driver_url('http://cida-test.er.usgs.gov/mda.lakes/')})
+
 clusterExport(c1, 'driver_fun')
 
-out = clusterApplyLB(c1, to_run, run_cal, years=1980:1999, driver_function=function(site_id){driver_fun(site_id, 'ACCESS')})
+out = clusterApplyLB(c1, to_run, run_cal, years=1980:1999, driver_function=driver_fun, datasource=NA)
+
 
 sprintf('%i lakes ran\n', sum(unlist(lapply(out, inherits, what='data.frame'))))
 
+
 ##results
 out_df = do.call('rbind', out[unlist(lapply(out, inherits, what='data.frame'))])
-write.table(out_df, 'c:/WiLMA/results/2016-03-19_ACCESS.tsv', sep='\t', row.names=FALSE)
+write.table(out_df, paste0('c:/WiLMA/results/', run_name, '.tsv'), sep='\t', row.names=FALSE)
+
+rmarkdown::render('demo/document_calibration_overview_template.Rmd', output_file=paste0(run_name,'.pdf'), 
+                  output_dir='c:/WiLMA/results/', params=list(out_df=out_df, run_message=run_comment))
+
 
 sqrt(mean((out_df$Observed_temp - out_df$Modeled_temp)^2, na.rm=TRUE))
 mean(out_df$Observed_temp - out_df$Modeled_temp, na.rm=TRUE)
 
-rmarkdown::render('demo/document_calibration_overview_template.Rmd', output_file='2016-03-19_ACCESS.pdf', 
-                  output_dir='c:/WiLMA/results/', params=list(out_df=out_df))
-
 
 # groups = split(to_run, ceiling(seq_along(to_run)/100))
 # out = list()

demo/document_calibration_overview_template.Rmd

@@ -5,10 +5,11 @@ date: "March 22, 2016"
 output: pdf_document
 params:
   out_df: !r data.frame()
+  run_message: !r "default message"
 
 ---
 
-Lets see how the run went. 
+Run description: `r params$run_message`
 
 
 
@@ -17,6 +18,8 @@ library(lubridate)
 library(glmtools)
 library(mda.lakes)
 
+out_df = params$out_df
+
 #out_df = read.table('c:/WiLMA/results/2016-03-21_NLDAS.tsv', header=TRUE, sep='\t', as.is=TRUE)
 Sys.setenv(tz='GMT')