-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathACM_GPP_ET_validation.r
231 lines (202 loc) · 11.3 KB
/
ACM_GPP_ET_validation.r
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
###
## function to return x,y coordinate from an array which is nearest to a provided lat / long value
###
closest2d <- function (id,lat,long,lat_in,long_in,nos_dim) {
# extract needed lat / long
lat1=lat_in[id] ; long1=long_in[id]
# calculate the distance between two points by Spherical law of the cosine
# mean radius of earth in km
R=6371 # 6378137 m (R source equitorial)
# convert degrees to radians
deg_to_rad=pi/180
# check lat long system
if (length(which(as.vector(long) > 180)) > 1) {stop("Input error closest2d: longitude should be -180 to +180")}
if (nos_dim == 1) {
## lat long are in single vectors repeating i.e. lat[1:10]=89.9,89.9,89.9... ; long[1:10]=-180,-160,-140....
# loop through to find the smallest distance
d_old=1e6
# check for locations which exactly coincide, complete calculatuion and finally remove NaN generated by value "1"
d = acos(sin(lat1*deg_to_rad)*sin(lat*deg_to_rad)+cos(lat1*deg_to_rad)*cos(lat*deg_to_rad)*cos((long*deg_to_rad)-(long1*deg_to_rad)))*R
match = which(is.na(d) == TRUE) ; d[match] = 0 #; print(match)
d = pmax(1e-6,d,na.rm=TRUE)
output=which(d == min(d))[1] ; d_old=d[output]
#if (d_old > 10) {print(paste("id = ",id," Minimum distance found (",round(d_old,digits=1),") is greater than 10 km from actual location"))}
rm(d_old)
} else if (nos_dim == 2) {
## lat and long are in two - dimensional arrays which co-varying
# loop through to find the smallest distance
d = sin(lat1*deg_to_rad)*sin(as.vector(lat)*deg_to_rad)+cos(lat1*deg_to_rad)*cos(as.vector(lat)*deg_to_rad)*cos((as.vector(long)*deg_to_rad)-(long1*deg_to_rad))
# check for locations which exactly coincide, complete calculatuion and finally remove NaN generated by value "1"
d = acos(d)*R ; match = which(is.na(d) == TRUE) ; d[match] = 0 #; print(match)
#if (min(d) > 10) {print(paste("id = ",id," Minimum distance found (",round(min(d),digits=1),") is greater than 10 km from actual location"))}
i=ceiling(which(d == min(d))/dim(lat)[1])
j=which(d == min(d))-floor(which(d == min(d))/dim(lat)[1])*dim(lat)[1]
output = list(j[1],i[1]) ; rm(i,j)
} else if (nos_dim == 3) {
## lat and long are in two 1-dimensional vectors but co-varying as in cartesian co-ordinates
# loop through to find the smallest distance
d_old=1e6
for (j in seq(1, length(lat))) {
# convert all to radians
d = acos(sin(lat1*deg_to_rad)*sin(lat[j]*deg_to_rad)+cos(lat1*deg_to_rad)*cos(lat[j]*deg_to_rad)*cos((long*deg_to_rad)-(long1*deg_to_rad)))*R
if (min(d) < d_old) {
possible_i = which(d == min(d))
#print(possible_i) ; print(j) ; print(length(long)) ; print(length(lat))
if (length(possible_i) > 1 & possible_i[1] == 1 & possible_i[length(possible_i)] == length(long)) {
possible_i = possible_i[which(long[possible_i] == long1)]
}
output = list(possible_i[1],j) ; d_old=min(d)
}
} # j loop
#if (d_old > 10) {print(paste("id = ",id," Minimum distance found (",round(d_old,digits=1),") is greater than 10 km from actual location"))}
rm(i,d_old)
}
# clean up
rm(d,R,lat1,long1) # ; gc() ; gc()
# return to the user
return(output)
} # end of function
###
## Create needed ACM_GPP_ET shared object
# Packages
library(RColorBrewer)
# set to the working directory that this script should be called from
setwd("/home/lsmallma/WORK/GREENHOUSE/models/ACM_GPP_ET/") ; wkdir = getwd()
# compile the shared object containing ACM_GPP and ACM_ET
system("gfortran ./src/ACM_GPP_ET.f90 ./src/ACM_GPP_ET_R_interface.f90 -o ./src/acm_gpp_et.so -fPIC -shared")
system("mv ./src/acm_gpp_et.so .")
###
## Borrow met data from an existing CARDAMOM analysis
drivers = read.csv("/home/lsmallma/gcel/ACM_GPP_ET_RECALIBRATION/output_files/acm_recal_with_spa_200pixels_continuous_timeseries_obs_iWUE_trunk_nowater_copy.csv")
###
## Define our output variables based on the grid of the CARDAMOM analysis we are borrowing
mean_lai = array(NA, dim=c(dim(drivers)[1]))
mean_gpp = array(NA, dim=c(dim(drivers)[1]))
mean_transpiration = array(NA, dim=c(dim(drivers)[1]))
mean_wetcanopyevap = array(NA, dim=c(dim(drivers)[1]))
mean_soilevaporation = array(NA, dim=c(dim(drivers)[1]))
mean_rootwatermm = array(NA, dim=c(dim(drivers)[1]))
mean_runoffmm = array(NA, dim=c(dim(drivers)[1]))
mean_drainagemm = array(NA, dim=c(dim(drivers)[1]))
mean_WUE = array(NA, dim=c(dim(drivers)[1]))
mean_wSWP = array(NA, dim=c(dim(drivers)[1]))
###
## Some ACM_GPP_ET parameters
output_dim=11 ; nofluxes = 8 ; nopools = 1 ; nopars = 4 ; nos_iter = 1
#soils_data=read.csv("/home/lsmallma/gcel/HWSD/processed_file/HWSD_sand_silt_clay_orgfrac_vector_with_lat_long_200.csv",header=TRUE)
# iterative process through the years...
for (n in seq(1,dim(drivers)[1])) {
if (n%%1000 == 0){print(paste("...beginning site:",n," of ",dim(drivers)[1], sep=""))}
# met note that the dimension here are different to that of drivers$met
met=array(-9999,dim=c(1,12))
met[,1] = 1 # day of analysis
met[,2] = drivers$sat_min[n] # min temperature (oC)
met[,3] = drivers$sat_max[n] # max temperature (oC)
met[,4] = drivers$swrad_avg[n]*(24*60*60*1e-6) # SW Radiation (MJ.m-2.day-1)
met[,5] = drivers$co2_avg[n] # CO2 ppm
met[,6] = drivers$doy[n] # day of year
met[,7] = drivers$ppt_avg[n]/(60*60) # rainfall (kgH2O.m-2.hr-1 -> kg.m-2.s-1)
met[,8] = -9999 # drivers$sat_avg[n] # avg temperature (oC)
met[,9] = drivers$wind_avg[n] # avg wind speed (m.s-1)
met[,10]= drivers$vpd_avg[n]*1000 # avg VPD (kPa->Pa)
# ecosystem state drivers now rather than meteorology
met[,11]= drivers$lai[n] # LAI
met[,12]= drivers$lai[n]*80 #100 # root C stocks
# parameters
parameters = array(NA, dim=c(nopars,nos_iter))
parameters[1,] = drivers$avgN[n] # foliar N (gN.m-2)
parameters[2,] = -9999 # min leaf water potential (MPa)
parameters[3,] = 100 # root biomass needed to reach 50 % depth
parameters[4,] = 2 # max root depth (m)
# other inputs
lat = drivers$lat[n]
# search location of soils data
# i1=unlist(closest2d(1,soils_data$lat_wanted,soils_data$long_wanted,drivers$lat[n],drivers$long[n],1))[1]
# soil_info=c(pmax(1,soils_data$sand_top[i1]),pmax(1,soils_data$sand_bot[i1]),pmax(1,soils_data$clay_top[i1]),pmax(1,soils_data$clay_bot[i1]) )
# if (soil_info[2] == 1) {soil_info[2] = soil_info[1]}
# if (soil_info[4] == 1) {soil_info[4] = soil_info[3]}
soil_info=c(drivers$sand_top[n],drivers$sand_bot[n],drivers$clay_top[n],drivers$clay_bot[n])
if (is.loaded("racmgppet") == FALSE) { dyn.load("./acm_gpp_et.so") }
tmp=.Fortran("racmgppet",output_dim=as.integer(output_dim),met=as.double(t(met)),pars=as.double(parameters)
,out_var=as.double(array(0,dim=c(nos_iter,(dim(met)[1]),output_dim)))
,lat=as.double(lat),nopars=as.integer(nopars),nomet=as.integer(dim(met)[2])
,nofluxes=as.integer(nofluxes),nopools=as.integer(nopools)
,nodays=as.integer(dim(met)[1])
,deltat=as.double(array(0,dim=c(as.integer(dim(met)[1])))),nos_iter=as.integer(nos_iter)
,soil_frac_clay=as.double(array(c(soil_info[3],soil_info[3],soil_info[4],soil_info[4]),dim=c(4)))
,soil_frac_sand=as.double(array(c(soil_info[1],soil_info[1],soil_info[2],soil_info[2]),dim=c(4))) )
output=tmp$out_var ; output=array(output, dim=c(nos_iter,(dim(met)[1]),output_dim))
if (n == dim(drivers)[1]) {dyn.unload("./acm_gpp_et.so")}
rm(tmp) ; gc()
# assign outputs to out final grids
mean_lai[n] = mean(output[,,1]) # lai
mean_gpp[n] = mean(output[,,2]) # GPP (gC.m-2.day-1)
mean_transpiration[n] = mean(output[,,3]) # transpiration (kg.m-2.day-1)
mean_wetcanopyevap[n] = mean(output[,,4]) # wetcanopy evaporation (kg.m-2.day-1)
mean_soilevaporation[n] = mean(output[,,5]) # soil evaporation (kg.m-2.day-1)
mean_wSWP[n] = mean(output[,,6]) # weighted soil water potential (MPa)
mean_WUE[n] = mean_gpp[n]/mean_transpiration[n] # water use efficiency (gC/kgH2O)
mean_rootwatermm[n] = mean(output[,,7]) # water in rooting zone (mm)
mean_runoffmm[n] = mean(output[,,8]) # surface runoff (mm)
mean_drainagemm[n] = mean(output[,,9]) # drainage from soil column (mm)
} # site loop
###
## Calculate some statistics
###
# R2
gpp_r2 = summary(lm(drivers$GPP~mean_gpp))$adj.r.squared
transpiration_r2 = summary(lm((drivers$Evap-drivers$soilevap-drivers$wetevap)~mean_transpiration))$adj.r.squared
soilevaporation_r2 = summary(lm(drivers$soilevap~mean_soilevaporation))$adj.r.squared
wetcanopyevap_r2 = summary(lm(drivers$wetevap~mean_wetcanopyevap))$adj.r.squared
# bias
gpp_bias = mean(mean_gpp-drivers$GPP, na.rm=TRUE)
transpiration_bias = mean(mean_transpiration-(drivers$Evap-drivers$soilevap-drivers$wetevap), na.rm=TRUE)
soilevaporation_bias = mean(mean_soilevaporation-drivers$soilevap, na.rm=TRUE)
wetcanopyevap_bias = mean(mean_wetcanopyevap-drivers$wetevap, na.rm=TRUE)
# rmse
gpp_rmse = sqrt(mean((drivers$GPP-mean_gpp)**2))
transpiration_rmse = sqrt(mean(((drivers$Evap-drivers$soilevap-drivers$wetevap)-mean_transpiration)**2))
soilevaporation_rmse = sqrt(mean((drivers$soilevap-mean_soilevaporation)**2))
wetcanopyevap_rmse = sqrt(mean((drivers$wetevap-mean_wetcanopyevap)**2))
print("GPP gC/m2/day")
print(gpp_r2) ; print(gpp_bias) ; print(gpp_rmse)
print("Transpiration kgH2O/m2/day")
print(transpiration_r2) ; print(transpiration_bias) ; print(transpiration_rmse)
print("Soil evaporation kgH2O/m2/day")
print(soilevaporation_r2) ; print(soilevaporation_bias) ; print(soilevaporation_rmse)
print("Wet canopy evaporation kgH2O/m2/day")
print(wetcanopyevap_r2) ; print(wetcanopyevap_bias) ; print(wetcanopyevap_rmse)
###
## Begin output to file
###
units=c("mean_lai = m2/m2","mean_gpp = gC/m2/day"
,"mean_transpiration = kgH2O/m2/day","mean_wetcanopyevap = kgH2O/m2/day","mean_soilevaporation = kg/H2O/m2/day"
,"mean_wSWP = MPa","mean_WUE = gC/kgH2O","mean_rootwatermm = kgH2O/m2"
,"mean_runoffmm = kgH2O/m2/day","mean_drainagemm = kgH2O/m2/day")
# Save output for later use
calibration_output = list(drivers = drivers,
units = units,
gpp_r2 = gpp_r2,
gpp_bias = gpp_bias,
gpp_rmse = gpp_rmse,
transpiration_r2 = transpiration_r2,
transpiration_bias = transpiration_bias,
transpiration_rmse = transpiration_rmse,
soilevaporation_r2 = soilevaporation_r2,
soilevaporation_bias = soilevaporation_bias,
soilevaporation_rmse = soilevaporation_rmse,
wetcanopyevap_r2 = wetcanopyevap_r2,
wetcanopyevap_bias = wetcanopyevap_bias,
wetcanopyevap_rmse = wetcanopyevap_rmse,
mean_lai = mean_lai,
mean_gpp = mean_gpp,
mean_transpiration = mean_transpiration,
mean_wetcanopyevap = mean_wetcanopyevap,
mean_soilevaporation = mean_soilevaporation,
mean_wSWP = mean_wSWP,
mean_WUE = mean_WUE,
mean_rootwatermm = mean_rootwatermm,
mean_runoffmm = mean_runoffmm,
mean_drainagemm = mean_drainagemm)
# Now save the file
save(calibration_output, file="./outputs/calibration_output.RData")