This document demonstrates the main functionality of the dWaSSIC package.
# function for checking libs f_lib_check<-function(libs){ for (lib in libs ){ if(lib %in% rownames(installed.packages())){ }else{ install.packages(lib,repos='http://cran.us.r-project.org') } } a<-lapply(libs, require, character.only = TRUE) } # load devtools package for installing github packages if("devtools" %in% installed.packages()){ library(devtools) }else{ install.packages("devtools") library(devtools) } # Install dWaSSIC package from github if("dWaSSI" %in% installed.packages()){ library(dWaSSIC) }else{ install_github("ln1267/dWaSSIC") install_github("tanerumit/sacsmaR") library(dWaSSIC) } # Load other required packages f_lib_check(c("raster","reshape2","parallel","lubridate","rgdal","rgeos","dplyr","ggplot2","leaflet","sacsmaR"))
Climate<-read.csv("INPUTS/Inputs_catchment/CLIMATE.TXT") LAI<-read.csv("INPUTS/Inputs_catchment/LANDLAI.TXT") SOIL<-read.csv("INPUTS/Inputs_catchment/SOILINFO.TXT") HUC<-read.csv("INPUTS/Inputs_catchment/CELLINFO.TXT") names(Climate)<-c("BasinID","Year","Month","Ppt_mm","Tavg_C") names(SOIL)[1:2]<-c("ID","BasinID") names(LAI)[1:2]<-c("BasinID","Year") names(HUC)[1:2]<-c("ID","BasinID") huc_lc_ratio<-as.matrix(HUC[c(4:length(HUC[1,]))]) hru_info<-HUC[c(1:3)] names(hru_info)<-c("BasinID","HRU_Lat","HRU_Lon") # Filter data based on climate data BasinIDs<-unique(Climate$BasinID) NBasins<-length(BasinIDs) Climate<-Climate[Climate$BasinID %in% BasinIDs,] HUC<-HUC[HUC$BasinID %in% BasinIDs,] LAI<-LAI[LAI$BasinID %in% BasinIDs,] SOIL<-SOIL[SOIL$BasinID %in% BasinIDs,]
# Set period for simulation Sim_dates<-list() Sim_dates[["Start"]] <- as.Date("2004/01/01") Sim_dates[["End"]] <- as.Date("2012/12/01" ) # Climate time range Sim_dates[["Start_climate"]] <- as.Date(paste0(min(Climate$Year),"/01/01")) Sim_dates[["End_climate"]] <- as.Date(paste0(max(Climate$Year),"/12/01")) # LAI time range Sim_dates[["Start_lai"]] <- as.Date(paste0(min(LAI$Year),"/01/01") ) Sim_dates[["End_lai"]] <- as.Date(paste0(max(LAI$Year),"/12/01")) Sim_dates[["Seq_date"]]<-seq.Date(as.Date(Sim_dates[["Start"]]),as.Date(Sim_dates[["End"]]),by = "month") Sim_dates[["Seq_date_climate"]]<-seq.Date(Sim_dates[["Start_climate"]],Sim_dates[["End_climate"]],by = "month") Sim_dates[["Sim_ind"]] <-which(Sim_dates[["Seq_date_climate"]] %in% Sim_dates[["Seq_date"]])
## Soil parameters par_sacsma = as.matrix(SOIL[,3:13] ) colnames(par_sacsma)<-toupper(colnames(par_sacsma)) ## PET parameters par_petHamon = rep(1,NBasins) ## Input climate data clim_prcp<-matrix(Climate$Ppt_mm,ncol=NBasins) clim_tavg<-matrix(Climate$Tavg_C,ncol=NBasins) ## Input LAI and ratio for each vegetation type huc_lc_ratio<-as.matrix(HUC[c(6:length(HUC[1,]))]) hru_lc_lai<-lapply(BasinIDs, function(x) as.matrix(subset(LAI,BasinID==x)[c(4:length(LAI[1,]))])) ## Hru info like Lat, long and elevation for each HUC hru_info<-HUC[c(1:5)] names(hru_info)<-c("ID","BasinID","HRU_Area","HRU_Lat","HRU_Lon")#,"HRU_Elev(m)","HRU_FlowLen(m)") hru_info["HRU_Elev(m)"]<-1000 hru_info["HRU_FlowLen(m)"]<-1000 # sample routing par par_routing<-stcroix$hru.par[1:22,28:31] # Load coefficients for SUN's ET and carbon calculation ET_coefs<-read.csv("INPUTS/Inputs_catchment/ET_theory.csv",nrows = length(huc_lc_ratio[1,])) names(ET_coefs)<-c("LC_ID","Intercept","P_coef","PET_coef","LAI_coef","P_PET_coef","P_LAI_coef","PET_LAI_coef","IGBP","LC_Name") WUE_coefs<-read.csv("INPUTS/Inputs_catchment/WUE_theory.csv",nrows = length(huc_lc_ratio[1,])) names(WUE_coefs)<-c("LC_ID","WUE","RECO_Interc", "RECO_Slope","IGBP","LC_Name")
flowR <- dWaSSIC(sim.dates=Sim_dates, warmup = 3,mcores = 1, par.sacsma = par_sacsma,par.petHamon=par_petHamon,par.routing =par_routing, hru.info = hru_info, clim.prcp = clim_prcp, clim.tavg = clim_tavg, hru.lai=NULL,hru.lc.lai=hru_lc_lai,huc.lc.ratio=huc_lc_ratio,WUE.coefs = WUE_coefs,ET.coefs = ET_coefs) out_weight<-lapply(c(1:length(flowR$HUC)),function (x) flowR$HUC[[x]]*HUC$Area_m2[x]/sum(HUC$Area_m2)) flowR_Catchment<-as.data.frame(sapply(names(out_weight[[1]]), function(var) apply(sapply(out_weight, function(x) x[[var]]),1,sum))) Rain_mean<-apply(clim_prcp, 1, mean,na.rm=T) df <- data.frame(date = Sim_dates$Seq_date,Srufflow = flowR$FLOW_SURF,Baseflow=flowR$FLOW_BASE,AET=flowR_Catchment$totaet,Rain=Rain_mean[Sim_dates$Sim_ind],nt=flowR_Catchment$tot) ggplot(df, aes(date, Srufflow)) + geom_line() + labs(x = "", y = "Surface flow (mm/month)") ggplot(df, aes(date, Baseflow)) + geom_line() + labs(x = "", y = "Base flow (mm/month)") ggplot(df, aes(date, Baseflow+Srufflow)) + geom_line() + labs(x = "", y = "Total flow (mm/month)") ggplot(df, aes(date, nt)) + geom_line() + labs(x = "", y = "Actual ET (mm/month)") ggplot(df, aes(date, Rain)) + geom_line() + labs(x = "", y = "Rainfall (mm/month)")
In this step the required original input includes:
1. Basin or subbasin boundaries, which can be created using data(DEM)
and ArcSwat toolbox in ArcGIS or using other watershed delineary tools.
2.
Basins<-readOGR("data/shps/Watersheds.shp") Basins$BasinID<-c(1:length(Basins[,1])) # get the contral coordinates of each polygon Basin_coords<-gCentroid(Basins, byid=TRUE) rownames(Basin_coords@coords)<-Basins$BasinID Basin_coords<-as.data.frame(Basin_coords) # Add latitude and longitude infor to the Basin if(!"lat" %in% names(Basins)) { Basins[["Lat"]]=Basin_coords$y Basins[["Long"]]=Basin_coords$x } # plot the basin leaflet(Basins) %>%addTiles() %>% addPolygons(color = "black", weight = 1, smoothFactor = 0.5)
LUCC_catchment<-hru_lc_ratio(classname ="data/Landcover/LUCC_Sun_IGBP.nc",shp = Basins,field = "BasinID") hru_lcs<-dcast(LUCC_catchment[c("BasinID","Class","Ratio")],BasinID~Class) hru_lcs[is.na(hru_lcs)]<-0.0 cellinfo<-Basins@data[c("BasinID","Shape_Area","Lat","Long","Elev","Len1")] cellinfo<-arrange(cellinfo,BasinID) cellinfo[c("Shape_Area","Elev","Len1")]<-round(cellinfo[c("Shape_Area","Elev","Len1")],0) cellinfo[c( "Lat","Long_")]<-round(cellinfo[c( "Lat","Long_")],3) cellinfo<-merge(cellinfo,hru_lcs,by="BasinID") cellinfo<-cbind("ID"=c(1:length(Basins)),cellinfo) names(cellinfo)[1:7]<-c("ID","BasinID","Area_m2","Lat","Long","Elev","Flowlen") write.table(cellinfo,"INPUTS/test/CELLINFO.TXT",sep = ',',row.names = FALSE)
Pre_catchment<-f_sta_shp_nc("data/Climate/Pre_MJ_mon_00_15.tif",Basins,varname = "P",start = 1982,zonal_field = "BasinID") Tmean_catchment<-f_sta_shp_nc("data/Climate/Temp_MJ_mon_00_15.tif",Basins,varname = "T",start = 1982,zonal_field = "BasinID") climate<-Pre_catchment climate$Tavg<-Tmean_catchment$T climate<-arrange(climate,BasinID,Year,Month) climate<-climate[c("BasinID","Year","Month","P","Tavg")] climate[c("P","Tavg")]<-round(climate[c("P","Tavg")],2) climate$BasinID<-as.integer(as.character(climate$BasinID)) names(climate)<-c("BasinID","Year","Month","Ppt_mm","Tavg_C") write.table(climate,"INPUTS/test/CLIMATE.TXT",row.names = F,sep=",")
climate_terr<-read.csv("data/Climate/MJ_subbasins_climate.csv") climate_terr$date<-as.Date(as.character(climate_terr$date),format="%Y%m%d") climate_terr[c("aet","pet","tmmn","soil","tmmx")]<-climate_terr[c("aet","pet","tmmn","soil","tmmx")]*0.1 climate_terr[c("pdsi")]<-climate_terr[c("pdsi")]*0.01 climate_terr$tmean<-(climate_terr$tmmn+climate_terr$tmmn)/2 climate_terr<-climate_terr[,-c(1,13)] names(climate_terr)[1]<-c("BasinID") climate_terr$Year<-as.integer(format(climate_terr$date,"%Y")) climate_terr$Month<-as.integer(format(climate_terr$date,"%m")) climate_terr<-arrange(climate_terr,BasinID,Year,Month) climate_terr[,c("pr","tmean")]<-round(climate_terr[,c("pr","tmean")],2) climate<-climate_terr[c("BasinID","Year","Month","pr","tmean")] names(climate)<-c("BasinID","Year","Month","Ppt_mm","Tavg_C") write.table(climate,"INPUTS/test/CLIMATE.TXT",row.names = F,sep=",") climate_terr%>% filter(date>="2000-01-01")%>% filter(BasinID %in% c("10","11"))%>% ggplot(aes(x=date,y=tmean))+geom_point()+geom_line()+facet_grid(BasinID~.)
f_sta_HRS_CMIP5<-function(i,shp){ tifs<-dir("/Dataset/backup/Climate_data/GMIP5/TIF/",pattern = "",full.names = T) tif_names<-dir("/Dataset/backup/Climate_data/GMIP5/TIF/",pattern = "") name<-substr(tif_names[i],1,10) da_stack<-raster(tifs[i]) print(range(da_stack)) #extract all data based on each shape boundary a<-raster::extract(da_stack,shp,fun=mean,df=T) rm(da_stack) c<-as.data.frame(t(a[-1])) # get the mean value for each shp feature names(c)<-shp@data$BasinID_Nu c } sta_CMIP5<-f_sta_HRS_CMIP5(1,Basins) f_Parallel_set("r") sta_CMIP5<-parLapply(cl,c(1:4536),f_sta_HRS_CMIP5,shp=HRS_basin) save(sta_CMIP5,file="data/HRS/sta_CMIP5.RData") sta_CMIP5<-do.call(rbind,sta_CMIP5) stopCluster(cl)
# this function is used for zonal LAI of each lc in the HRU hru_lai<-hru_lc_zonal(classname = "data/Landcover/LUCC_Sun_IGBP.nc", daname = "data/LAI/LAI_1982_2013.tif", shp = Basins, field = "BasinID") lcs<-paste0("Class_",names(hru_lcs)[-1]) f_fillmatix<-function(a,lcs){ a<-round(a,2) prel<-length(a[1,])+1 if(prel< length(lcs)+1){ lacks<-lcs[which(! lcs %in% colnames(a))] for (i in c(1:length(lacks))) a<-cbind(a,lcadd=0) colnames(a)[prel:length(a[1,])]<-lacks a<-a[,lcs] } } ha<-lapply(hru_lai, f_fillmatix,lcs) hru_lais<-do.call(rbind,ha) hru_lais<-cbind("ID"=rep(as.integer(names(hru_lai)),each=length(hru_lai[[1]][,1])), "Year"=rep(c(1982:2014),each=12), "Month"=c(1:12), hru_lais) hru_lais<-as.data.frame(hru_lais) hru_lais<-arrange(hru_lais,ID,Year,Month) write.table(hru_lais,"INPUTS/test/LANDLAI.TXT",sep = ',',row.names = FALSE)
SOIL<-brick("data/Soil/SOIL_BNU.nc") SOIL_catchment<-NA SOIL_catchment<-extract(SOIL,Basins,fun=mean,na.rm=T,df=T) SOIL_catchment<-SOIL_catchment[-1] names(SOIL_catchment)<-c("uztwm", "uzfwm" , "uzk", "zperc" , "rexp" , "lztwm" , "lzfsm", "lzfpm", "lzsk" , "lzpk" , "pfree") # fill NA values for (i in c(1:length(SOIL_catchment))) SOIL_catchment[[i]][is.infinite(SOIL_catchment[[i]])]<-NA for (i in c(1:length(SOIL_catchment))) SOIL_catchment[[i]][is.na(SOIL_catchment[[i]])]<-30 SOIL_catchment<-cbind(BasinID=Basins$BasinID,SOIL_catchment) SOIL_catchment[c(2:12)]<-round(SOIL_catchment[c(2:12)],2) write.table(SOIL_catchment[c(1:12)],"INPUTS/test/SOILINFO.TXT",sep = ',',row.names = FALSE)
generalinfo<-c() generalinfo[1]<-"Minjiang watershed Simulation by Cell or Catchment, Sep. 2018 (Ning Liu)" generalinfo[2]<-"0 Scenarios (See Explaination Below)" generalinfo[3]<-"0 Catchment or Grid scale (0 for Catchment, 1 for Grid)" generalinfo[4]<-"22, 62, 1951, 10 (No watersheds, NoYEARS, first year of climate data, number of land cover categories)" generalinfo[5]<-"2000, 2012 (Start and end years for LAI data)" generalinfo[6]<-"5,2000, 2012 (Years for warmup, and Start and end years for simulation)" generalinfo[7]<-"0.0 (FOREST REDUCTION TO CONVERT CROPLAND (FRACTION 0-1))" generalinfo[8]<-"0.0 (Reduction fraction of Leaf Area Index, 0.0=no change, 1.0= convert to bared soil) " generalinfo[9]<-"0.0 INTIAL SNOWPACK (MM)" generalinfo[10]<-"" generalinfo[11]<-"" generalinfo[12]<-"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" generalinfo[13]<-"Simulation Scenarios" generalinfo[14]<-"0 - Historic baseline 1901-2002" generalinfo[15]<-"1 - Deforestation" generalinfo[16]<-"2 - LAI reduction (same land use, but different cover) " writeLines(generalinfo,"INPUTS/test/GENERAL.TXT")
water<-read.csv("~/Downloads/FOR_NING/out/CONUSF2F_MA_WATER.TXT") routpar<-read.csv("~/Downloads/FOR_NING/INPUTS/HUC12_CONUS_ROUTE_TABLE.txt") water$flow<-water$YLD_ALL*water$AREA_M2/1000/1000000 ha<-hrurouting(water,routpar,mc_cores = 5) out<-read.csv("~/Downloads/FOR_NING/out/CONUSF2F_MA_FLOW_HUC12.TXT") head(out[,c(1,5)],10) head(ha[,c(1,15)],10)
cd src/ make rm ../OUTPUTS/* ./dWaSSIC 1 ../INPUTS/Inputs_catchment ../OUTPUTS #./dWaSSIC 1 ../INPUTS/Inputs_grid ../OUTPUTS echo "finished" make clean
# Input data U <- 2 #m/s L <- 25000 #m x <- 0.2 delta_t <- 6*60*60 #s outflow.init <- 12 #m^3/s muskingum <- function(inflow, U, L, x, delta_t, outflow.init){ # Muskingum method of flood routing. Requires inflow hydrograph, returns # outflow hydrograph. Assumes kt = L/U. kt <- L/U c0 <- (-kt*x + 0.5*delta_t)/(kt - kt*x + 0.5*delta_t) c1 <- (kt*x + 0.5*delta_t)/(kt - kt*x + 0.5*delta_t) c2 <- (kt - kt*x - 0.5*delta_t)/(kt - kt*x + 0.5*delta_t) outflow <- rep(0, length(inflow)) outflow[1] <- outflow.init for (i in 2:length(inflow)){ outflow[i] <- c0*inflow[i] + c1*inflow[i-1] + c2*outflow[i-1] } return(outflow) } df$date<-as.Date(df$date) df["flow"] <- muskingum(df$Srufflow, U, L, x, delta_t, outflow.init) ggplot(data= df,aes(x=df$date, y=df$flow)) + geom_line() + geom_point() + scale_x_continuous("Time (days)") + scale_y_continuous("Discharge (cubic meters per second)") + ggtitle("Flood Hydrograph 25km Downstream") + theme_bw()
opt_ns<-function(a,obs,sim){ abs(sum(a*obs-sim)) } xx=c(1,2,3,4,5,6,7,8,9,10) zz=xx*5 optimize(f = opt_ns,c(0,10),tol = 0.0001,obs=xx, sim=zz) require(graphics) f <- function (x, a) (x - a)^2 xmin <- optimize(f, c(0, 1), tol = 0.0001, a = 1/3) xmin ## See where the function is evaluated: optimize(function(x) x^2*(print(x)-1), lower = 0, upper = 10) ## "wrong" solution with unlucky interval and piecewise constant f(): f <- function(x) ifelse(x > -1, ifelse(x < 4, exp(-1/abs(x - 1)), 10), 10) fp <- function(x) { print(x); f(x) } plot(f, -2,5, ylim = 0:1, col = 2) optimize(fp, c(-4, 20)) # doesn't see the minimum optimize(fp, c(-7, 20)) # ok optim(par=c(0), fn=opt_ns, obs=xx, sim=yy) yy=c(30,40,22,33,40) funk=function(param,x,y,z){ a=rep(param[1],5) b=param[2] d=param[3] fit=sum((y-(a+b*x+z*d))^2) return(fit) } optim(par=c(1,1,1), fn=funk, x=xx, y=yy, z=zz)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.