monitoring/R_code_faPAR.r at master · mooradios/monitoring · GitHub

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
# How to look at a chemical cycling from a satellite

# Set working directory
# setwd("~/lab/")
# setwd("/Users/utente/lab") #mac
setwd("C:/lab/") # windows

# require libraries

library(raster)
library(rasterVis)
library(rasterdiv)

plot(copNDVI)

# Reclassify the copNDVI for removing the blue part (water)
copNDVI <- reclassify(copNDVI, cbind(253:255, NA))

plot(copNDVI)

levelplot(copNDVI)

# import the data (aggregated already)
faPAR10 <- raster("faPAR10.tif")

levelplot(faPAR10)


#export PDF

pdf("copNDVI.pdf")
levelplot(copNDVI)
dev.off()

pdf("faPAR10.pdf")
levelplot(faPAR10)
dev.off()

#############################

#2nd day

# Set working directory
# setwd("~/lab/")
# setwd("/Users/utente/lab") #mac
setwd("C:/lab/") # windows

# load the workspace

load(".RData")

ls()

faPAR10

#require libraries

library(raster)

library(rasterdiv)

library(rasterVis)

# the original file of fapar from Copernicus is 2 GB

# let's see how much space is needed for 8-BIT set


# see the range of file copNDVI

copNDVI

# write a raster file

writeRaster(copNDVI, "copNDVI.tif") #5.3 MB, that's why we are using this files instead of original radiance files

levelplot(faPAR10)

######################
#3rd day

# regression model between faPAR and NDVI

# create 2 variables
erosion <- c(12, 14, 16, 24, 26, 40, 55, 67) # kg per square meter
hm <- c(30, 100, 150, 200, 260, 340, 460, 600) # ppm

plot(erosion, hm, col="red", pch=19, xlab ="erosion", ylab = "heavy metals", cex=2)


# see the relation by using the linear model

model1 <- lm(hm ~ erosion)


summary(model1) # to check the r squared, p-value and etc.


# include the linear model into the plot

abline(model1)

# Set working directory
# setwd("~/lab/")
# setwd("/Users/utente/lab") #mac
setwd("C:/lab/") # windows

# load the workspace

load(".RData")


#require libraries

library(raster)

library(rasterdiv)

library(rasterVis)

library(sf)


# plot the faPAR10

plot(faPAR10)

# plot the copNDVI

plot(copNDVI)

# check the number of cells in the faPAR10
faPAR10

# random points from the image

random.points <- function(x,n)
{
  lin <- rasterToContour(is.na(x))
  pol <- as(st_union(st_polygonize(st_as_sf(lin))), 'Spatial') # st_union to dissolve geometries
  pts <- spsample(pol[1,], n, type = 'random')
}

# generate random points from the image

pts <- random.points(faPAR10, 1000)

plot(faPAR10)

points(pts, col="red", pch=19)

# passing values from the map to the sample

copNDVIp <- extract(copNDVI, pts)
faPAR10p <- extract(faPAR10, pts)

copNDVIp # to see the 1000 random points

faPAR10p

dev.off()

# build the model
# photosynthesis vs biomass

model2 <- lm(faPAR10p ~ copNDVIp)

# summary of the model

summary(model2)

plot(copNDVIp, faPAR10p, col="green", xlab="biomass", ylab="photosynthesis")
abline(model2, col="red")
# plot shows the previous patterns that we observed (check the levelplots of copNDVI and faPAR for more information)

levelplot(copNDVI)
levelplot(faPAR10)