# 06/19/2023 # run ech2o on mtbs fires, starting with rice ridge # modify ech2o spatial function for this task # code here is from ech2o_web_points_v0.12_multipoint WRITE_HOURLY_NCDF=FALSE WRITE_TIFF=TRUE # define spatial res of all outputs. This is important... spatial_res=240 RUN.SPATIAL=TRUE RUN.CLIMATE=TRUE RUN.BINARY.CLIMATE=TRUE OVERWRITE.CLIMATE=FALSE BUILD.CONFIG=TRUE RUN.ECH20=FALSE # Add run type here: pre fire, post-fire, uniform. This gets passed to spatial_function. RUN_TYPE="PREFIRE" # PREFIRE, POSTFIRE are other options require(terra) require(sf) require(foreach) require(doParallel) # Jan 24th another check of missing stations. base.dir="/mnt/DataDrive1/data/ech2o_projects/ech2o_ET/" ncdf.dir <- "/mnt/DataDrive1/data/airtemps/ncdf/" prcp.dir <- "/mnt/DataDrive1/data/PRISM_precip/daily_prcp/" srad.dir <- "/mnt/mordor3/data/airtemps/ncdf/" soil.dir <- "/mnt/DataDrive1/data/soil_grids/adjusted_grids_v2/" aux.dir <- "/mnt/DataDrive1/data/ech2o_projects/points_westUS_2023/aux_files/" nc.dir="/mnt/DataDrive1/data/ech2o_projects/ech2o_mtbs/ech2o_ncdf_240/" # use RAP forest cover vcf.dir="/mnt/DataDrive1/data/RAP_cover/" # source some basic functions (listfiles, fdates) and ech2o support functions. source("/mnt/DataDrive1/data/ech2o_projects/ech2o_mtbs/scripts/my functions.R") #source("/mnt/DataDrive2/data/watershed_modeling/scripts/echo_spatial_function_30m.R") source("/mnt/DataDrive1/data/ech2o_projects/ech2o_mtbs/scripts/echo_spatial_function_user_DEM_06192023.R") source("/mnt/DataDrive1/data/ech2o_projects/ech2o_mtbs/scripts/ech2o_daily_to_hourly_functions.R") mtpoly <- st_read("/mnt/DataDrive1/data/MTBS_zh/mtbs_perimeter_2026/bsp_perims_DD.shp") if(spatial_res==240) { echo.dir=paste0(base.dir, "ech2o_outputs_240/") ech2o.ncdf=paste0(base.dir, "ech2o_ncdf_240/") } if(spatial_res==100) { echo.dir=paste0(base.dir, "ech2o_outputs_100/") ech2o.ncdf=paste0(base.dir, "ech2o_ncdf_100/") } run_type="HUC" REPORT_SOIL=1; REPORT_TSKIN=1; REPORT_FLOW=0; REPORT_ET=1 # try using a square of data for fire mask.watershed="TRUE" # define run interval; always 3-hourly ObsPerDay <- 8 # define the hours for which temperature will be retrieved do_hours <- seq.int(0,23, by=24/ObsPerDay) INIT_YEARS=4; POST_YEARS=5 def250 <- rast("/mnt/DataDrive1/data/waterbalance_250/def_1981-2010_annual_mean.tif") #require(future) #future::plan(callr, workers=16) #cl=makeCluster(12) #registerDoParallel(cl) #iFire=438 # Rice Ridge #for(iFire in 1:length(id.list)) { #fire_id=id.list[iFire] #foreach(fire_id=id.list) %dopar% { library(sf) library(terra) library(ncdf4) # must install package from website (http://cirrus.ucsd.edu/~pierce/ncdf/) def250 <- rast("/mnt/DataDrive1/data/waterbalance_250/def_1981-2010_annual_mean.tif") # Lolo Peak fire fire_id="MT4667411426820170715" # Rice Ridge #fire_id="mt4726811348520170724" PID=fire_id out.dir=paste0(echo.dir, fire_id, "/") dir.create(out.dir) # create raster mask from mtbs polygon testID <- which(mtpoly$event_id==fire_id)[1] poly1=mtpoly[testID,] fire_date=as.character(poly1$ig_date) fire_year=as.numeric(substr(fire_date, 1,4)) dnbr_tmp=vect(poly1) rast_tmplt=rast(dnbr_tmp, res=c(.00027, .00027), crs='epsg:4326'); rast_tmplt=setValues(rast_tmplt, 1:ncell(rast_tmplt)) dnbr=rasterize(dnbr_tmp, rast_tmplt) dnbr_proj=project(dnbr, 'epsg:4326', method="near") writeRaster(dnbr_proj, file=paste0(out.dir, "/", "polyrast_", fire_id, ".tif"), overwrite=T) fire_ext=ext(dnbr_proj) fire_ext=extend(fire_ext, c(2,2)) #dnbr=crop(dnbr, fire_ext) dnbr_wgs=project(dnbr_proj, 'epsg:4326') ext_wgs <- ext(dnbr_wgs) huc1=dnbr_wgs # using this old object to avoid replacing it many times below latitude=ext_wgs[3] # define start/end dates with init and post year vars above start_year=fire_year-INIT_YEARS ; end_year=fire_year+POST_YEARS if(start_year < 1986) { start_year <- 1986 } if(end_year > 2024) { end_year <- 2024 } start.date <- paste0(start_year, "-10-01"); end.date <- paste0(end_year, "-09-30") do.dates <- seq(as.Date(start.date), as.Date(end.date), "days") # covert dates to Year/month for targeting netcdf monthly files. yearmonth <- unique(format(do.dates, "%Y%m")) # create the date fields used for hourly to daily conversion here. sdate <- as.POSIXct(paste0(start.date, " 00:00"), format="%Y-%m-%d %H:%M") edate <- as.POSIXct(paste0(end.date, " 21:59"), format="%Y-%m-%d %H:%M") dt <- seq(sdate, edate, by="3 hours") date <- substr(dt, 1,10) spatial.dir <- paste(out.dir, "/Spatial/", sep="") if(!exists(spatial.dir)) {dir.create(spatial.dir)} climate.dir <- paste(out.dir, "/Climate/", sep="") if(!exists(climate.dir)) {dir.create(climate.dir)} if(RUN_TYPE=="PREFIRE") { result.dir <- paste(out.dir, "/Results/", sep="") } if(RUN_TYPE=="POSTFIRE") { result.dir <- paste(out.dir, "/Results_post/", sep="") } if(RUN_TYPE=="UNIFORM") { result.dir <- paste(out.dir, "/Results_uni/", sep="") } if(!exists(result.dir)) {dir.create(result.dir)} # call spatial function here if(RUN.SPATIAL==TRUE) { # check spat files here #splist=length(list.files(spatial.dir, pattern=".map")) #if(splist<50){ fulldem <- rast("/mnt/DataDrive1/data/landfire/us_dem.tif") #dnbr=project(dnbr, 'epsg:5070', res=c(240,240), method="near") ech2o_spatial_function_huc(PID, dnbr, spatial_res) #} } ################################## # climate input section project.dir <- out.dir if(RUN.CLIMATE==TRUE) { demc <- rast(paste0(spatial.dir, PID, "_dem_NOBUFFER.tif")) NAflag(demc) <- -9999 # demcp=project(demc, crs(dnbr_wgs)) # dnbr_mask=which(dnbr==1) #dnbr_mask=terra::project(dnbr, 'epsg:5070', res=c(240,240), method="near") dnbr_mask=terra::project(dnbr, demc, method="near", align=F) demc=mask(demc, dnbr_mask, maskvalues=NA) #dem.mask <- Which(!is.na(demc)) # make slightly larger extent object for cropping crop_extent <- ext_wgs crop_extent[3] <- crop_extent[3]-.01; crop_extent[4] <- crop_extent[4] +.01; crop_extent[1] <- crop_extent[1]-.01; crop_extent[2] <- crop_extent[2] +.01; crop_extent_prcp=extend(crop_extent, .05) # t1 <- Sys.time() # print("starting Tmin extraction") # # ncdf crop method (should be faster) # ncfile1 <- paste0(ncdf.dir, "conus_tmin_", "200101", "_finescale.nc") # ncb1 <- rast(ncfile1) # tmin.test <- crop(ncb1, crop_extent) # tmin.rast <- subset(tmin.test, 1) #nc1 <- nc_open(ncfile1) #lonidx <- which(nc1$dim$lon$vals>xmin(huc1) & nc1$dim$lon$vals<= xmax(huc1)) #latidx <- which(nc1$dim$lat$vals>ymin(huc1) & nc1$dim$lat$vals<= ymax(huc1)) #nc_close(nc1) #tmin.array <- c() tmin.stack <- rast() for(iFile in 1:length(yearmonth)) { ncfile1 <- paste0(ncdf.dir, "conus_tmin_", yearmonth[iFile], "_finescale.nc") if(!file.exists(ncfile1)) { ncfile1 <- paste0(ncdf.dir, "conus_tmin_", yearmonth[iFile], ".nc") } #nc1 <- nc_open(ncfile1) #mydat.all <- ncvar_get(nc1, varid= 'tmin', start=c(lonidx[1], latidx[1],1), count=c(ncol(tmin.test),nrow(tmin.test), -1)) nc1=rast(ncfile1) nc1c=terra::crop(nc1, crop_extent) rm(nc1) #nc_close(nc1) #tmin.array <- abind(tmin.array, mydat.all) tmin.stack=c(tmin.stack, nc1c) rm(nc1c) } #tmin.stack <- brick(nrows=nrow(tmin.rast), ncols=ncol(tmin.rast), xmn=xmin(tmin.rast), xmx=xmax(tmin.rast), ymn=ymin(tmin.rast), ymx=ymax(tmin.rast), nl=dim(tmin.array)[3], crs=projection(tmin.rast)) #tmin.stack <- setValues(tmin.stack, tmin.array) #tmin.stack <- rast(tmin.array, extent=ext(tmin.rast), crs=crs(tmin.rast)) #t2 <- Sys.time() #tdiff <- t2-t1 #print("finished tmin extraction") #print(tdiff) #tmin.size <- object.size(tmin.array) #print(tmin.size) #rm(tmin.array) gc() print("starting Tmax extraction") tstart <- Sys.time() tmax.stack <- rast() for(iFile in 1:length(yearmonth)) { ncfile1 <- paste0(ncdf.dir, "conus_tmax_", yearmonth[iFile], ".nc") if(!file.exists(ncfile1)) { ncfile1 <- paste0(ncdf.dir, "conus_tmax_", yearmonth[iFile], "_finescale.nc") } nc1=rast(ncfile1) nc1c=terra::crop(nc1, crop_extent) rm(nc1); #nc_close(nc1) #tmin.array <- abind(tmin.array, mydat.all) tmax.stack=c(tmax.stack, nc1c) } # extract minimum RH print("starting RH extraction") tstart <- Sys.time() #rmin.array <- c() rmin.stack <- rast() for(iFile in 1:length(yearmonth)) { ncfile1 <- paste0(ncdf.dir, "conus_rhmin_", yearmonth[iFile], "_8s.nc") if(!file.exists(ncfile1)) { ncfile1 <- paste0(ncdf.dir, "conus_rhmin_", yearmonth[iFile], "_finescale.nc") } if(!file.exists(ncfile1)) { ncfile1 <- paste0(ncdf.dir, "conus_rhmin_", yearmonth[iFile], ".nc") } nc1=rast(ncfile1) nc1c=terra::crop(nc1, crop_extent) rm(nc1) #nc_close(nc1) #tmin.array <- abind(tmin.array, mydat.all) rmin.stack=c(rmin.stack, nc1c) } #rmin.stack <- rast(rmin.array, extent=ext(tmin.rast), crs=crs(tmin.rast)) rmin.stack <- rmin.stack/100 # t2 <- Sys.time() # tdiff <- t2-tstart # print("finished rhmin extraction") # print(tdiff) #rm(rmin.array) gc() # rhmax print("beginning RHmax extraction") tstart=Sys.time() rmax.stack <- rast() for(iFile in 1:length(yearmonth)) { ncfile1 <- paste0(ncdf.dir, "conus_rhmax_", yearmonth[iFile], "_8s.nc") if(!file.exists(ncfile1)) { ncfile1 <- paste0(ncdf.dir, "conus_rhmax_", yearmonth[iFile], "_finescale.nc") } if(!file.exists(ncfile1)) { ncfile1 <- paste0(ncdf.dir, "conus_rhmax_", yearmonth[iFile], ".nc") } nc1=rast(ncfile1) nc1c=terra::crop(nc1, crop_extent) rm(nc1) #nc_close(nc1) #tmin.array <- abind(tmin.array, mydat.all) rmax.stack=c(rmax.stack, nc1c) } #rmax.stack <- rast(rmax.array, extent=ext(tmin.rast), crs=crs(tmin.rast)) rmax.stack <- rmax.stack/100 testmin=rmin.stack <=.01 rmin.stack[testmin]<- .01 testmin=rmax.stack >=1 rmax.stack[testmin] <- 1 t2 <- Sys.time() tdiff <- t2-tstart print("finished rhmax extraction") print(tdiff) #rm(rmax.array) gc() # prism precip. Will need to resample here. prcp.list <- list.files(prcp.dir, pattern=".nc", full=T) start.yr <- as.numeric(substr(start.date, 1,4))-1 end.yr <- as.numeric(substr(end.date, 1,4)) do.years <- seq(start.yr, end.yr, 1) prcp.dates <- seq(as.Date(paste0(start.yr, "-01-01")), as.Date(paste0(end.yr, "-12-31")), "days") prcp.stack <- rast() for(iYear in 1:length(do.years)) { ncfile1 <- paste0(prcp.dir, "prism_precipitation_4km_", do.years[iYear], ".nc") nc1 <- rast(ncfile1) prcp.year <- crop(nc1, crop_extent_prcp) rm(nc1) prcp.stack <- c(prcp.stack, prcp.year) } # retrieve just the dates we need get.start <- match(start.date, as.character(prcp.dates)) get.stop <- match(end.date, as.character(prcp.dates)) prcp.stack <- subset(prcp.stack, get.start:get.stop) # # Solar radiation #srad.array <- c() srad.stack=rast() for(iFile in 1:length(yearmonth)) { ncfile1 <- paste0(srad.dir, "conus_srad_", yearmonth[iFile], ".nc") if(!file.exists(ncfile1)) { ncfile1 <- paste0(srad.dir, "conus_srad_", yearmonth[iFile], "_finescale.nc") } nc1=rast(ncfile1) nc1c=terra::crop(nc1, crop_extent) srad.stack=c(srad.stack, nc1c) rm(nc1) } #rm(srad.array) gc() print("finished Temperature RH SRAD and Precip extraction for HUC") ######################## # write the stacked geotiffs out here? Probably need to convert to ncdf # to retain the date fields? # do this if we want to provide way to extract gridded weather spatially. # after writing, mask all grids and project them. # Oct 16th 2023. write hourly ncdf on demand. for nelson calcs ################################## # begin daily to hourly conversion for HUC run here # Start with Tmin/Tmax # system("ping -w 10 1.1.1.1") print(dim(demc)) print(dim(tmin.stack)) print("non-NA DEM cells") print(length(which(is.na(values(demc))))) #print("saving workspace for testing") #save.image(paste0(final.dir, "workspace.RData")) print("projecting weather raster stacks") print(Sys.time()) tmin.stack <- terra::project(tmin.stack, demc) print("finished projecting tmin weather rasters") tmax.stack <- terra::project(tmax.stack, demc) print("finished projecting tmax weather rasters") rmin.stack <- terra::project(rmin.stack, demc) print("finished projecting rmin weather rasters") rmax.stack <- terra::project(rmax.stack, demc) print("finished projecting rmax weather rasters") prcp.stack <- terra::project(prcp.stack, demc) print("finished projecting prcp weather rasters") srad.stack <- terra::project(srad.stack, demc) print("finished projecting srad weather rasters") print("finished projecting all weather rasters") print(Sys.time()) # create wind rasters (all 1 for now) print("creating wind raster stack") wind.stack <- setValues(tmin.stack, 1) # make sure there are no extra NA's in precip grids. if(mask.watershed) { print("masking weather raster stacks") tmin.stack <- mask(tmin.stack, demc, maskvalues=NA) tmax.stack <- mask(tmax.stack, demc, maskvalues=NA) rmin.stack <- mask(rmin.stack, demc, maskvalues=NA) rmax.stack <- mask(rmax.stack, demc, maskvalues=NA) prcp.stack <- mask(prcp.stack, demc, maskvalues=NA) srad.stack <- mask(srad.stack, demc, maskvalues=NA) wind.stack <- mask(wind.stack, demc, maskvalues=NA) } if(WRITE_TIFF==TRUE) { tiff.dates=strftime(do.dates, format="%Y%m%d") names(tmin.stack) <- tiff.dates names(tmax.stack) <- tiff.dates names(rmax.stack) <- tiff.dates names(rmin.stack) <- tiff.dates names(srad.stack) <- tiff.dates names(prcp.stack) <- tiff.dates tiff.year=substr(tiff.dates,1,4) keep.tiff.yr=which(tiff.year==fire_year) tmin.stack.tiff=subset(tmin.stack, keep.tiff.yr) tmax.stack.tiff=subset(tmax.stack, keep.tiff.yr) srad.stack.tiff=subset(srad.stack, keep.tiff.yr) rmin.stack.tiff=subset(rmin.stack, keep.tiff.yr) rmax.stack.tiff=subset(rmax.stack, keep.tiff.yr) prcp.stack.tiff=subset(prcp.stack, keep.tiff.yr) writeRaster(tmin.stack.tiff, file=paste0(out.dir, "tmin_", fire_year, ".tif"), overwrite=T) writeRaster(tmax.stack.tiff, file=paste0(out.dir, "tmax_", fire_year, ".tif"), overwrite=T) writeRaster(rmin.stack.tiff, file=paste0(out.dir, "rmin_", fire_year, ".tif"), overwrite=T) writeRaster(rmax.stack.tiff, file=paste0(out.dir, "rmax_", fire_year, ".tif"), overwrite=T) writeRaster(srad.stack.tiff, file=paste0(out.dir, "srad_", fire_year, ".tif"), overwrite=T) writeRaster(prcp.stack.tiff, file=paste0(out.dir, "prcp_", fire_year, ".tif"), overwrite=T) } # create longwave stack here. # try with arrays. 2.7 minutes for 1 year. Much faster than raster::stack but still slow #tmin.a <- as.array(tmin.stack); tmax.a <- as.array(tmax.stack); rmin.a <- as.array(rmin.stack); rmax.a <- as.array(rmax.stack) # test here with a single vector of values from each brick. This takes less than 10 seconds! print("creating downward longwave raster stack") print(Sys.time()) # can this be done in raster form? #tavg=((tmin.stack+273.15) + (tmax.stack+273.15))/2 #rhavg=(rmin.stack+rmax.stack)/2 #dlr.stack=dlr_fun_clear(tavg, rhavg) tmin.val <- values(tmin.stack); tmax.val <- values(tmax.stack); rmin.val <- values(rmin.stack); rmax.val <- values(rmax.stack) temp.val <- (tmin.val+273.15 + tmax.val+273.15)/2; rh.val <- (rmin.val+rmax.val)/2 dlr.val <- dlr_fun_clear(temp.val, rh.val) dlr.stack <- setValues(tmin.stack, dlr.val) #print(Sys.time()) ############################# if(RUN.BINARY.CLIMATE==TRUE) { print("converting daily data to hourly for HUC") all_hours <- seq.int(0,23, by=24/ObsPerDay) # create binary header file info here vals <- values(demc) vals <- vals[which(!is.na(vals))] vals <- round(vals, digits=1) rnum <- length(do.dates)*ObsPerDay # how many rows in climate input file: number of time periods for tmin rnames <- seq(1, rnum, 1) cnum <- length(vals) cnames <- seq(1, length(vals), 1) # cnames <- format(cnames, scientific=F, quote=F) # proper format for each header line of binary file rnum <- as.integer(rnum) rnames <- as.numeric(rnames) cnum <- as.integer(cnum) head.info <- vector(mode="character", length=256) # write the first 4 lines of file, then close # working with hourly, Tmin, Max and Avg should be redundant. maybe write one file, then just copy it? inputfile.nameM <- paste(climate.dir, "tmin", ".bin", sep="") inputfile.nameX <- paste(climate.dir, "tmax", ".bin", sep="") inputfile.nameA <- paste(climate.dir, "tavg", ".bin", sep="") temp.ncdf=rast() zzM <- file(inputfile.nameM, "wb") #open an output file connection writeBin(head.info, zzM, size=1) # exactly 256 character cnames <- as.integer(cnames) writeBin(rnum, zzM, size=4) writeBin(rnames, zzM, size=4) writeBin(cnum, zzM, size=4) writeBin(cnames, zzM, size=4) # begin daily to hourly binary. # look at using arrays or vectors; raster is slow. # no projection or masking should be needed. # begin loop through days here, converting to hourly and writing hourly data to binary file. i=1 # define the hours for which temperature will be retrieved #for(i in 1:2) { for(i in 1:length(do.dates)) { if(i==1) { do_hours <- seq.int(0,24, by=24/ObsPerDay)[(length(all_hours)/2+2):length(all_hours)] } else { do_hours <- all_hours } #print(i) # print(paste(Sys.time(), "temperature day ", i)) JDay <- strptime(do.dates[i], format="%Y-%m-%d")$yday rastm <- subset(tmin.stack, i) rastx <- subset(tmax.stack, i) if(i==1) { prevmin <- rastm; prevmax <- rastx } else { prevmin <- subset(tmin.stack, i-1) prevmax <- subset(tmax.stack, i-1) } if(i < length(do.dates)) { nextmin <- subset(tmin.stack, i+1) } else { nextmin <- rastm } Tmin <- values(rastm); Tmax <- values(rastx); PrevMin <- values(prevmin); PrevMax <- values(prevmax); NextMin <- values(nextmin) # loop through hours in day, estimate temperature at hour iHr for(iHr in 1:length(do_hours)) { hourcount <- do_hours[iHr] hourtemp <- temperature_hourly(hourcount, latitude, JDay, Tmin, Tmax, PrevMin, PrevMax, NextMin) if(WRITE_HOURLY_NCDF==TRUE) { r1=setValues(rastm, hourtemp) temp.ncdf=c(temp.ncdf, r1) } hourtemp <- as.vector(na.omit(hourtemp)) #print(summary(hourtemp)) #print(Sys.time()) # hour_rast <- overlay(lat_rast, JDay_rast, prevmin, prevmax, rastm, rastx, nextmin, fun=my_daily_to_hourly_fun1 ) # hour_vals <- values(hour_rast) writeBin(hourtemp, zzM, size=4) } } # end of loop through days in current year close(zzM) # copy the hourly file created above and rename as Tavg and Tmax file.copy(inputfile.nameM, inputfile.nameX, overwrite=T) file.copy(inputfile.nameM, inputfile.nameA, overwrite=T) #close(zzX) #close(zzA) # end writing temperaure binary files hourly ######################################### # # End of temperature binary file creation ############################################# ################################################# ################################################ # Daily to hourly humidity for HUC run print("converting daily RH data to hourly for HUC") inputfile.nameRH <- paste(climate.dir, "rhavg", ".bin", sep="") # add hourly ncdf file here for RH rh.ncdf=rast() zzRH <- file(inputfile.nameRH, "wb") #open an output file connection writeBin(head.info, zzRH, size=1) # exactly 256 character writeBin(rnum, zzRH, size=4) writeBin(rnames, zzRH, size=4) writeBin(cnum, zzRH, size=4) writeBin(cnames, zzRH, size=4) # begin loop through days here, converting to hourly and writing hourly data to binary file. # note that we reverse RHmin and RHmax here so the hourly curves are correct. t1=Sys.time() for(i in 1:length(do.dates)) { JDay <- strptime(do.dates[i], format="%Y-%m-%d")$yday + 1 rastx <- subset(rmin.stack, i) rastm <- subset(rmax.stack, i) if(i==1) { prevmin <- rastm; prevmax <- rastx } else { prevmax <- subset(rmin.stack, i-1) prevmin <- subset(rmax.stack, i-1) } if(i < length(do.dates)) { nextmin <- subset(rmax.stack, i+1) } else { nextmin <- rastx } rmin <- values(rastm); rmax <- values(rastx); PrevMin <- values(prevmin); PrevMax <- values(prevmax); NextMin <- values(nextmin) # loop through hours in day, estimate temperature at hour iHr # add 4 hour shift forward here if it's the first day of run. if(i==1) { do_hours <- seq.int(0,24, by=24/ObsPerDay)[(length(all_hours)/2+2):length(all_hours)] } else { do_hours <- all_hours } for(iHr in 1:length(do_hours)) { hourcount <- do_hours[iHr] hourtemp <- humidity_hourly(hourcount, latitude, JDay, rmin, rmax, PrevMin, PrevMax, NextMin) hourtemp[hourtemp>1] <- 1 hourtemp[hourtemp<0] <- .01 if(WRITE_HOURLY_NCDF==TRUE) { r1=setValues(rastm, hourtemp) rh.ncdf=c(rh.ncdf, r1) } hourtemp <- as.vector(na.omit(hourtemp)) #print(paste0(hourcount, " ", hourtemp)) writeBin(hourtemp, zzRH, size=4) } } # end of loop through days in current year close(zzRH) # end daily to hourly RH ############## #t2=Sys.time(); print(t2-t1) ###################### # daily to hourly SRAD for HUC run print("converting daily srad to hourly for HUC") inputfile.nameRad <- paste(climate.dir, "swrad", ".bin", sep="") if(WRITE_HOURLY_NCDF==TRUE & RUN.CLIMATE==TRUE) { srad.ncdf=rast() } zzR <- file(inputfile.nameRad, "wb") #open an output file connection writeBin(head.info, zzR, size=1) # exactly 256 character cnames <- as.integer(cnames) writeBin(rnum, zzR, size=4) writeBin(rnames, zzR, size=4) writeBin(cnum, zzR, size=4) writeBin(cnames, zzR, size=4) for(i in 1:length(do.dates)) { JDay <- strptime(do.dates[i], format="%Y-%m-%d")$yday + 1 rad1 <- subset(srad.stack, i) netradval <- values(rad1) # loop through hours in day, estimate temperature at hour iHr if(i==1) { do_hours <- seq.int(0,24, by=24/ObsPerDay)[(length(all_hours)/2+2):length(all_hours)] hours.day1=length(do_hours) } else { do_hours <- all_hours } for(iHr in 1:length(do_hours)) { hourcount <- do_hours[iHr] hourrad <- srad_hourly_huc(hourcount, latitude, JDay, netradval) if(WRITE_HOURLY_NCDF==TRUE) { r1=setValues(rad1, hourrad) srad.ncdf=c(srad.ncdf, r1) } hourrad <- as.vector(na.omit(hourrad)) writeBin(hourrad, zzR, size=4) } # end loop through hours in day } # end loop through days in year close(zzR) # End daily to hourly radiation for HUC #################### ##################### # begin precip daily to hourly conversion # write the first 4 lines of file, then close inputfile.nameP <- paste(climate.dir, "precip", ".bin", sep="") prcp.ncdf=rast() zzP <- file(inputfile.nameP, "wb") #open an output file connection writeBin(head.info, zzP, size=1) # exactly 256 character writeBin(rnum, zzP, size=4) writeBin(rnames, zzP, size=4) writeBin(cnum, zzP, size=4) writeBin(cnames, zzP, size=4) #for(i in 1:2) { for(iDay in 1:length(do.dates)) { #print(iDay) if(iDay==1) { do_hours <- seq.int(0,24, by=24/ObsPerDay)[(length(all_hours)/2+2):length(all_hours)] } else { do_hours <- all_hours } # JDay <- strptime(dates_year[i], format="%Y%m%d")$yday prcp1 <- subset(prcp.stack, iDay) #testNA <- Which(is.na(prcp1)) #prcp1[testNA] <- 0 #prcp1 <- mask(prcp1, dem.mask, maskvalue=0) prcp.vals.full <- (values(prcp1))/(1000*86400) # convert to meters/sec prcp.vals <- as.vector(na.omit(prcp.vals.full)) for(iHr in 1:length(do_hours)) { if(WRITE_HOURLY_NCDF==TRUE) { r1=setValues(prcp1, prcp.vals.full*1000) prcp.ncdf=c(prcp.ncdf, r1) } writeBin(prcp.vals, zzP, size=4) } } # end loop through days in year close(zzP) # End precip file write for hourly data ############################################### ################################################## # downward longwave radiation inputfile.nameLW <- paste(climate.dir, "lwrad", ".bin", sep="") zzLW <- file(inputfile.nameLW, "wb") #open an output file connection writeBin(head.info, zzLW, size=1) # exactly 256 character writeBin(rnum, zzLW, size=4) writeBin(rnames, zzLW, size=4) writeBin(cnum, zzLW, size=4) writeBin(cnames, zzLW, size=4) for(i in 1:length(do.dates)) { vals <- as.vector(na.omit(values(subset(dlr.stack, i)))) if(i==1) { do_hours <- seq.int(0,24, by=24/ObsPerDay)[(length(all_hours)/2+2):length(all_hours)] } else { do_hours <- all_hours } for(iHr in 1:length(do_hours)) { writeBin(vals, zzLW, size=4) } } close(zzLW) # end of downward longwave binary file creation for HUC run ######################################### ####################################### # wind speed inputfile.nameW <- paste(climate.dir, "wind", ".bin", sep="") zzW <- file(inputfile.nameW, "wb") #open an output file connection writeBin(head.info, zzW, size=1) # exactly 256 character writeBin(rnum, zzW, size=4) writeBin(rnames, zzW, size=4) writeBin(cnum, zzW, size=4) writeBin(cnames, zzW, size=4) for(i in 1:length(do.dates)) { # here define # of days to process by length ofSWR. avoid leap year probs. vals <- as.vector(na.omit(values(subset(wind.stack, i)))) if(i==1) { do_hours <- seq.int(0,24, by=24/ObsPerDay)[(length(all_hours)/2+2):length(all_hours)] } else { do_hours <- all_hours } for(iHr in 1:length(do_hours)) { writeBin(vals, zzW, size=4) } } close(zzW) # end wind daily to hourly ################## print("done converting daily data to hourly for HUC") } # end RUN.BINARY ################## # end climate file generation for HUC run } # end RUN.CLIMATE ################################ # } # end OVERWRITE.CLIMATE #################### #################### # write out hourly NCDF files here for fuel moisture calcs # prepare and write multiband, multilayer ncdf. # if(WRITE_HOURLY_NCDF & RUN.CLIMATE==TRUE) { # time(srad.ncdf) <- dt[6:length(dt)] # time(temp.ncdf) <- dt[6:length(dt)] # time(rh.ncdf) <- dt[6:length(dt)] # time(prcp.ncdf) <- dt[6:length(dt)] # fm.ncdf=sds(temp.ncdf, rh.ncdf, srad.ncdf, prcp.ncdf) # names(fm.ncdf) <- c("tmp", "rh", "srad", "prcp") # writeCDF(fm.ncdf, file=paste0(climate.dir, "fm_3hourly.nc") ) #} ########################## # build config file and run ech2o here ############################ if(BUILD.CONFIG=="TRUE") { print("building ech2o config file") setwd(out.dir) config.file <- "config_ech2o_web_gridrun.ini" config.base <- readLines(paste(aux.dir, config.file, sep="")) config.outfile <- paste("config_", PID, ".ini", sep="") ndays=length(do.dates) # NOTE: To get max. tskin returned as daily output, we are eliminating # the 0,3, and 6 hours of an 8x daily run. This makes the last grid of each day the 12 pm hour. # make sure the simulation hours and # of output daily grids align with date sequence! time0_line <- paste("Simul_end = ", ndays*86400-5*(24/ObsPerDay*3600) ) time1_line <- paste("Simul_tstep = ", 24/ObsPerDay*3600, sep="") time2_line <- paste("Clim_input_tstep = ", 24/ObsPerDay*3600, sep="") #time3_line <- paste("Report_interval = ", 24/ObsPerDay*3600, sep="") #time3_line <- paste("Report_interval = ", 86400, sep="") time3_line <- paste("Report_interval = ", 10800, sep="") config.base[43] <- time0_line config.base[44] <- time1_line config.base[45] <- time2_line config.base[46] <- time3_line tmin.line <- paste("MinAirTemp = ", "tmin", ".bin", sep="") tmax.line <- paste("MaxAirTemp = ", "tmax", ".bin", sep="") tavg.line <- paste("AirTemperature = ", "tavg", ".bin", sep="") prcp.line <- paste("Precipitation = ", "precip", ".bin", sep="") rh.line <- paste("RelativeHumidity = ", "rhavg", ".bin", sep="") wind.line <- paste("WindSpeed = ", "wind", ".bin", sep="") swrad.line <- paste("IncomingShortWave = ", "swrad", ".bin", sep="") lwrad.line <- paste("IncomingLongWave = ", "lwrad", ".bin", sep="") config.base[61] <- prcp.line config.base[62] <- tavg.line config.base[63] <- tmax.line config.base[64] <- tmin.line config.base[65] <- rh.line config.base[66] <- wind.line config.base[67] <- lwrad.line config.base[68] <- swrad.line if(REPORT_FLOW==1 & run_type=="HUC"){ flow.line <- "Report_Streamflow = 1" config.base[148] <- flow.line } if(REPORT_SOIL==1 & run_type=="HUC"){ flow.line <- "Report_Soil_Water_Content_L1 = 1" config.base[150] <- flow.line } if(REPORT_TSKIN==1 & run_type=="HUC"){ flow.line <- "Report_Skin_Temperature = 1" config.base[162] <- flow.line } if(REPORT_ET==1 & run_type=="HUC"){ flow.line <- "Report_Transpiration = 1" config.base[179] <- flow.line } write.table(config.base, file=paste(out.dir, config.outfile, sep=""), col.names=F, row.names=F, quote=F, append=F, sep="\t ") ####################################### # end writing config file ####################################### rm(dlr.stack); rm(rmax.stack); rm(rmin.stack); rm(tmin.stack); rm(tmax.stack); rm(fulldem) gc() } # end BUILD.CONFIG #} # end loop through fires. #stopCluster(cl)