-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathEdHulPOL_development.txt
161 lines (139 loc) · 7.42 KB
/
EdHulPOL_development.txt
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
###
## Development for the EdHulPOL model
###
# The EdHulPOL model is a R based pollen dispersion model.
# The model allows for calibration of model parameters
# given known contempary data and subsequent application
# of these parameters for paleo-ecological studies
# to assess the senstivity of lakes to change in ecosystem composition.
# Novel in this application is the propogation of parameter uncertainties
# and iteration of the hypothetical land surfaces. Therefore, robust
# probability based estimates can be made (i.e. using Monte Carlo approaches).
# EdHulPOL is a modified version of the HULPOL v4 model.
# Bunting & Middleton (2005) Modelling pollen dispersal and deposition using HUMPOL software, including simulating windroses and irregular lakes,
# Review of Paleobotany & Palynology, vol 134, 184-196
## Model options
deposition_model=3 # select which deposition model?
## INPUTS
nos_species=2 # number of taxa considered
lake_size=1400 # radius of lake (m)
domain_radius=500*0.5 # radius of domain (pixels)
resolution=200 # resolution / pixel size (m)
fall_speed=c(0.18,0.250) # fall speed (m.s-1) of pollen; 0 = internal calculation
wind_sp=3 # wind speed (m.s-1)
radius=10e-6 # radius of pollen grain (m)
pollen_productivity=c(1,1.5) # runif(1,0,15)
## Parameters
gamma=0.125 # turbulence parameter = 2*gamma
diffusion_coefficient=0.12 # 0.3 # Vertical diffusion coefficient (m^1/8)
gravity=9.81 # gravity (m s-2)
rho0=2e6 # particle density (g m-3)
rho=1.27e6 # air density density (g m-3)
mu=1.8e-2 # dynamic viscosity (g m-1 s-1)
# species distance stuff using % covers
sp_hold=read.table("/home/lsmallma/WORK/EdHumPol/forest_example.txt")[,1]
sp_hold2=read.table("/home/lsmallma/WORK/EdHumPol/savanna_example.txt")[,1]
land_cover=read.table("/home/lsmallma/WORK/EdHumPol/smear_example2.txt",header=T)
water=land_cover$water
sp_hold=land_cover$forest
sp_hold2=land_cover$savanna
temp1=sp_hold/(sp_hold+sp_hold2+water)
temp2=sp_hold2/(sp_hold+sp_hold2+water)
sp_hold=temp1*100 ; sp_hold[which(is.na(sp_hold))]=0
sp_hold2=temp2*100 ; sp_hold2[which(is.na(sp_hold2))]=0
# what are my distance increments
distance_from_centre=seq(0,domain_radius*resolution,resolution) # basically radius of spatial domain and resolution (m)
distance_from_centre=distance_from_centre[-1]
# what is the total extend in steps
extent=(domain_radius*resolution*2)/resolution
# create vectorised versions of the 2D space
i_extent=rep(1:extent,times=extent) ; j_extent=rep(1:extent, each=extent)
# calculate distance from the centre of the grid
centre_i=median(1:extent) ; centre_j=median(1:extent)
distance_array=sqrt((abs(i_extent-centre_i)**2)+(abs(j_extent-centre_j)**2))*resolution
# which of these are wls()ithin the lake area, assuming circular lake
water_locations=which(distance_array < lake_size)
nos_water=length(water_locations)
##
### Calculate the mean species abundance at a given distance from the lake
# First allocate the distance variables and the nos_species dimension
mean_sp_abundance_at_distance=array(0, dim=c(length(distance_array),nos_species))
# next we loop through the input file to assess our distance vectorised
for (i in seq(1, length(sp_hold))) {
# periodic update
if (ceiling((i/length(sp_hold))*100) == (i/length(sp_hold))*100) {print(paste("Species maps ",round((i/length(sp_hold))*100,1)," complete %",sep=""))}
# actual calculation
mean_sp_abundance_at_distance[which(distance_array < (i*resolution) & distance_array > ((i-1)*resolution)),1] = sp_hold[i]
mean_sp_abundance_at_distance[which(distance_array < (i*resolution) & distance_array > ((i-1)*resolution)),2] = sp_hold2[i]
}
# update the distance vector
distance_array=array(distance_array, dim=c(length(distance_array),(nos_water)))
for (w in seq(1, nos_water)) {
# periodic updates
if (ceiling((w/nos_water)*100) == (w/nos_water)*100) { print(paste("Distance maps ",round((w/nos_water)*100,1)," % complete",sep="")) }
# locate new center, i.e. the water pixel
centre_i=water_locations[w]%%extent ; centre_j=ceiling(water_locations[w]/extent)
# now create multiple distance maps for each of the water locations
distance_array[,w]=sqrt((abs(i_extent-centre_i)**2)+(abs(j_extent-centre_j)**2))*resolution
}
# clean up
i_extent = 0 ; j_extent = 0 ; rm(i_extent,j_extent)
# ensure a minimum distance (1 m) in distance arrays
distance_array=as.vector(distance_array)
distance_array[which(distance_array == 0)] = 1
# generate spatial structures now
distance_array=array(distance_array, dim=c(extent,extent,nos_water))
mean_sp_abundance_at_distance=array(mean_sp_abundance_at_distance, dim=c(extent,extent,nos_species))
# check the result
par(mfrow=c(2,3)) ; library(fields)
image.plot(distance_array[,,1], main="Distance from centre (m)")
image.plot(mean_sp_abundance_at_distance[,,1], main="Forest species abundance (%)")
image.plot(mean_sp_abundance_at_distance[,,2], main="Savanna species abundance (%)")
# match matrix size with the species distribution information
pollen_deposition=array(NA, dim=c(dim(mean_sp_abundance_at_distance),nos_water))
if (deposition_model==1) {
for (w in seq(1, nos_water)) {
pollen_deposition[,,1:nos_species,w]=distance_array[,,w]**-1
}
} else if (deposition_model==2) {
pollen_deposition[,,1:nos_species,w]=distance_array[,,w]**-2
} else if (deposition_model==3) {
for (i in seq(1, nos_species)) {
print(paste("Species ",i," of ",nos_species,sep=""))
# Presentice model (Sutton 1953; Prentice 1988)
if (fall_speed[i] == 0) {fall_speed[i]=(2*(radius**2)*gravity*(rho0-rho))/(9*mu)}
beta_sp=(4*fall_speed[i])/((2*gamma)*wind_sp*(pi**0.5)*diffusion_coefficient)
# loop through the water tiles, possibly better to vectorise this also
for (w in seq(1, nos_water)) {
# periodic update
if (ceiling(((w*i)/(nos_water*nos_species))*100) == ((w*i)/(nos_water*nos_species))*100) {
print(paste("...water pixel = ",w," of ",nos_water,sep=""))
print(paste("...deposition matrix calculation ", round(((w*i)/(nos_water*nos_species))*100,1)," % complete",sep=""))
}
pollen_deposition[,,i,w]=beta_sp*gamma*(distance_array[,,w]**(gamma-1))*exp(-beta_sp*distance_array[,,w]**gamma)
} # water loop
} # species loop
} # all conditions
# clean up
distance_array = 0 ; rm(distance_array)
beta_sp = 0 ; rm(beta_sp)
# check the result
#image.plot(pollen_deposition[,,1,1], main="pollen deposition (Forest)")
#image.plot(pollen_deposition[,,2,1], main="pollen deposition (Savanna)")
# Prentice-Sugita model (Sugita 1994) modified by Bunting & Middleton (2005); where pollen load is linearly related to distance weighted plant abundance
pollen_load_at_lake_prediction=array(NA, dim=c(nos_species,nos_water))
for (w in seq(1,nos_water)) {
pollen_load_at_lake_prediction[1,w]=sum(pollen_productivity[1]*mean_sp_abundance_at_distance[,,1]*pollen_deposition[,,1,w])
pollen_load_at_lake_prediction[2,w]=sum(pollen_productivity[2]*mean_sp_abundance_at_distance[,,2]*pollen_deposition[,,2,w])
}
# number of pollen grains expected; mean across all water pixels
print(rowMeans(pollen_load_at_lake_prediction))
print(rowSums(pollen_load_at_lake_prediction))
plot(pollen_load_at_lake_prediction[1,], main="Pollen loads sp 1")
plot(pollen_load_at_lake_prediction[2,], main="Pollen loads sp 2")
# clean up again
pollen_load_at_lake_prediction = 0 ; rm(pollen_load_at_lake_prediction)
mean_sp_abundance_at_distance = 0 ; rm(mean_sp_abundance_at_distance)
pollen_deposition = 0 ; rm(pollen_deposition)
# garbage collection
gc(reset=TRUE, verbose=FALSE)