Introduction

This document demonstrates the main functionality of the dWaSSIC package.

load required document

# 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"))

Test run

Read input

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

# 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"]])

Convert input parameters

## 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")

Run the simulation using dWaSSIC funtion

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)")

Input data process

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.

Read Catchments boundary

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)

CELLINFO.TXT

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)

CLIMATE.TXT

Zonal from local data

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 data from GEE's TERRACLIMATE

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~.)

CMIP5 (optional)

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)

LANDLAI.TXT

# 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)

SOILINFO.TXT

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)

GENERAL.TXT

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")

Routing

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)

Fortran version

RUN fortran version

cd src/
make 
rm  ../OUTPUTS/*
./dWaSSIC  1 ../INPUTS/Inputs_catchment ../OUTPUTS
#./dWaSSIC  1 ../INPUTS/Inputs_grid ../OUTPUTS
echo "finished"
make clean

Underdevelopment

Testing muskingum Routing

# 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()

Testing optimise

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) 


ln1267/dWaSSI documentation built on Dec. 3, 2019, 4:39 a.m.