-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmetriche.Rmd
361 lines (302 loc) · 18.4 KB
/
metriche.Rmd
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
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
---
title: "metriche"
author: "Enrico Gabrielli F.A.R.M. Facilitazioni Agroecologiche Regionali Mobili"
date: "17/1/2023"
bibliography: vindicta_samc.bib
output:
md_document:
variant: markdown_github
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
# VINDICTA progetto 16.1 PSR Emilia-Romagna 2014-2020
## azione 3, analisi agroecosistema
### Calcolo delle metriche di analisi dell'agroecosistema
I grafici delle metriche sono dei plot a "ridge", ossia curve:
- la y, ovvero l'altezza delle curve, indica il numero di celle del raster
- x è il valore della metrica, per la resistenza in scala normale, per le altre metriche in scala logaritmica in base 10.
La scala dei colori utilizzata corrisponde a quella delle mappe.
```{r pacchetti, message=FALSE, include=FALSE}
library(terra) # per analizzare il raster
library(raster)
library(scales)
library(tidyverse)
library(ggplot2)
library(terrainr)
library(stringr)
library(dplyr)
library(listr)
library(ggridges)
library(patchwork)
library(gridExtra)
library(RColorBrewer)
```
```{r variabili, include=FALSE}
cartella <- './cartografia/output' #cartella in cui si trovano i raster in formato tif
```
```{r funzioni, include=FALSE}
fun_nomi_plot <- function(plot,nomi) {plot + ggtitle(nomi) +
theme(legend.position = "none",
legend.title = element_blank(), #rimuovi testi che non servono come "value" in legenda
legend.key.width = unit(1.5,'cm'), #change legend key width
plot.title = element_text(size=10), #size del titolo
axis.text = element_text(size = 7),
axis.text.x = element_text( angle = 90)
)} #funzione per aggiungere il titolo ai plot già creati
```
```{r raster files, echo=FALSE, message=FALSE, warning=FALSE}
raster_dir <- list.dirs(path = cartella)
raster_files <- list.files(path=cartella, #lista di file tif (nelle cartelle ci sono anche i README.md)
pattern = "\\.tif$",
full.names = TRUE,
recursive = T)
rasters <- lapply(raster_files, rast) #carico tutti i raster
rast_df <- lapply(rasters,function(x) terra::as.data.frame(x,xy=T)) #per poter fare ggplot trasformo tutto in data.frame
estrai <- strsplit(raster_files, "(\\\\|/)+") #si poteva fare anche con strextract in un solo passaggio, ma non ci sono riuscito
raster_files <- sapply(estrai, function(x) paste0(x[4],'-',x[5]))
for (i in seq_along(rast_df)){ #questo loop serve solo a rinominare i dataframe nella lista, importante perché poi i plot prenderanno qui il titolo
names(rast_df) <- raster_files
}
rast_df <- lapply(rast_df, setNames, c('x','y','value')) #la colonna dei dataframe deve essere uguale
#medie <- lapply(rast_df,function(x) mean(x[,3]))
rast_df_br <- bind_rows(rast_df,.id = "data_frame") #creo un'unica tabella con due colonne, una con il nome del raster, l'altra con i valori
rast_df_br$data_frame <- as.factor(rast_df_br$data_frame)
```
```{r foto_aeree, message=FALSE, warning=FALSE, include=FALSE}
cartella_foto_aeree <- './cartografia/fotointerpretazione/foto_aeree'
foto_aeree_files <- list.files(path=cartella_foto_aeree, #lista di file tif
pattern = "\\.tif$",
full.names = TRUE,
recursive = T)
wms_raster <- mapply(function(x) stack(x),foto_aeree_files) #carico la foto aerea nelle 3 componenti RGB - per una questione di semplicità non assegno l'estensione sul crs
fun_plot_rgb <- function(raster) ggplot() +
geom_spatial_rgb(
data = raster,
mapping = aes(
x=x,
y=y,
r=red,
g=green,
b=blue
)) +
labs(x = NULL, y = NULL)+ #tolgo i label agli assi
coord_fixed() #per avere un quadrato sempre
plot_foto_aeree <- lapply(wms_raster,fun_plot_rgb)
nomi_plot <- as.list(str_extract(names(wms_raster),"([^\\/]+)(?=\\.\\w+$)")) #per aggiungere i titoli per ogni plot faccio una lista dei nomi dei raster-dataframe
plot_foto_aeree <- mapply(plot=plot_foto_aeree,nomi=nomi_plot,fun_nomi_plot,SIMPLIFY = F) #ho due liste, una dei plot e una dei nomi, applico una funzione multipla, simplify=F serve a ricreare due liste, altrimenti per ogni layer dei plot usa la funzione
```
```{r heatmap resistenza, message=FALSE, include=FALSE}
resistenza_df <- rast_df[grepl("resistenza", names(rast_df))] #estraggo solo i dataframe della resistanza
fun_plot_resistance <- function(x) ggplot(x,aes(x=x,y=y))+ #funzione per l'heatmap
geom_tile(aes(fill=value))+
scale_fill_gradientn(colours = rev(brewer.pal(7, "RdYlGn")),
# trans='log',#scala non logaritmica per la resistenza
label = function(x) sprintf("%.2f", x),
limits=c(0.01, 1) #per avere le legende dei vari plot confrontabili
) +
theme(axis.text.x.bottom = element_blank(),
axis.ticks.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank(),
legend.position = "none",
legend.title = element_blank()
)+
labs(x = NULL, y = NULL) + #tolgo i label agli assi
coord_fixed() #per avere un quadrato sempre
plot_resistance <- lapply(resistenza_df,fun_plot_resistance)
#nomi_plot <- as.list(str_extract(names(resistenza_df),"([^\\-]+)(?=\\.\\w+$)")) #per aggiungere i titoli per ogni plot faccio una lista dei nomi dei raster-dataframe
#plot_resistance <- mapply(plot=plot_resistance,nomi=nomi_plot,fun_nomi_plot,SIMPLIFY = F) #ho due liste, una dei plot e una dei nomi, applico una funzione multipla, simplify=F serve a ricreare due liste, altrimenti per ogni layer dei plot usa la funzione
plot_resistance_dc <- wrap_plots(rev(plot_resistance),ncol = 1)
```
```{r plot resistenza, fig.ext='png', message=FALSE, warning=FALSE, include=FALSE}
#utilizziamo il ridge plot per una questione tecnica: con ggbeeswarm in questo caso non riesco a visualizzare una proporzione corretta tra i numeri di punti-pixel
#per la resistenza non uso la scala logaritmica
resistenza <- filter(rast_df_br,str_detect(data_frame,"resistenza")) #estraggo i dati di resistenza
resistenza$data_frame <- as.factor(str_extract(resistenza$data_frame,"([^\\-]+)(?=\\.\\w+$)")) #rinomino i fattori
ridge_resistance <- ggplot(resistenza,aes(x=value,y=data_frame,height = stat(density),fill=stat(x)))+
geom_density_ridges_gradient(stat = "density",rel_min_height = 0.01,scale=0.9)+
scale_fill_gradientn(colours = rev(brewer.pal(7, "RdYlGn")),
label = function(x) sprintf("%.2f", x),
limits=c(0.001, 1) #per avere le legende dei vari plot confrontabili
) +
scale_y_discrete(expand = c(0.01, 0)) +
labs(x=NULL,y = NULL) +
theme(legend.key.height = unit(50, 'points'),
legend.position = "left")
```
```{r heatmap assorbimento, message=FALSE, include=FALSE}
assorbimento_df <- rast_df[grepl("assorbimento", names(rast_df))] #estraggo solo i dataframe dell'assorbimento
fun_plot_absorbance <- function(x) ggplot(x,aes(x=x,y=y))+ #funzione per l'heatmap
geom_tile(aes(fill=value))+
scale_fill_gradientn(colours = rev(brewer.pal(7, "RdYlGn")),
trans='log10',
label = function(x) sprintf("%.2f", x)
,limits=c(0.000001, 1) #per avere le legende dei vari plot confrontabili. 10-5 è un valore empirico
) +
theme(axis.text.x.bottom = element_blank(),
axis.ticks.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank(),
legend.position = "none",
legend.title = element_blank()
)+
labs(x = NULL, y = NULL) + #tolgo i label agli assi
coord_fixed() #per avere un quadrato sempre
plot_absorbance <- lapply(assorbimento_df,fun_plot_absorbance)
#nomi_plot <- as.list(str_extract(names(assorbimento_df),"([^\\-]+)(?=\\.\\w+$)")) #per aggiungere i titoli per ogni plot faccio una lista dei nomi dei raster-dataframe
#plot_absorbance <- mapply(plot=plot_absorbance,nomi=nomi_plot,fun_nomi_plot,SIMPLIFY = F) #ho due liste, una dei plot e una dei nomi, applico una funzione multipla, simplify=F serve a ricreare due liste, altrimenti per ogni layer dei plot usa la funzione
plot_absorbance_dc <- wrap_plots(rev(plot_absorbance),ncol = 1)
```
```{r plot assorbanza, fig.ext='png', message=FALSE, warning=FALSE, include=FALSE}
assorbimento <- filter(rast_df_br,str_detect(data_frame,"assorbimento")) #estraggo i dati di assorbimento
assorbimento$data_frame <- as.factor(str_extract(assorbimento$data_frame,"([^\\-]+)(?=\\.\\w+$)")) #rinomino i fattori
ridge_assorbance <- ggplot(assorbimento,aes(x=value,y=data_frame,height = stat(count),fill=stat(x)))+
geom_density_ridges_gradient(stat = "density",rel_min_height = 0.01,scale=0.9)+
scale_fill_gradientn(colours = rev(brewer.pal(7, "RdYlGn")),
trans = "log10",
# label = function(x) sprintf("%.2f", x),
breaks=trans_breaks('log10', function(x) 10^x),
labels=trans_format('log10', math_format(10^.x))
,limits=c(0.000001, 1) #per avere le legende dei vari plot confrontabili
) +
scale_x_continuous(
trans='log10'
,breaks=trans_breaks('log10', function(x) 10^x),
labels=trans_format('log10', math_format(10^.x))
) +
scale_y_discrete(expand = c(0.01, 0)) +
annotation_logticks(sides = "b")+
labs(x=NULL,y = NULL) +
theme(legend.key.height = unit(50, 'points'),
legend.position = "left")
```
```{r heatmap dispersione, message=FALSE, include=FALSE}
dispersione_df <- rast_df[grepl("dispersione", names(rast_df))] #estraggo solo i dataframe della dispersione
fun_plot_dispersion <- function(x) ggplot(x,aes(x=x,y=y))+ #funzione per l'heatmap
geom_tile(aes(fill=value))+
scale_fill_gradientn(colours = brewer.pal(7, "RdYlGn"),
trans='log10',
label = function(x) sprintf("%.2f", x),
limits=c(1, 100) #per avere le legende dei vari plot confrontabili
) +
theme(axis.text.x.bottom = element_blank(),
axis.ticks.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank(),
legend.position = "none",
legend.title = element_blank()
)+
labs(x = NULL, y = NULL) + #tolgo i label agli assi
coord_fixed() #per avere un quadrato sempre
plot_dispersion <- lapply(dispersione_df,fun_plot_dispersion)
#nomi_plot <- as.list(str_extract(names(dispersione_df),"([^\\-]+)(?=\\.\\w+$)")) #per aggiungere i titoli per ogni plot faccio una lista dei nomi dei raster-dataframe
#plot_dispersion <- mapply(plot=plot_dispersion,nomi=nomi_plot,fun_nomi_plot,SIMPLIFY = F) #ho due liste, una dei plot e una dei nomi, applico una funzione multipla, simplify=F serve a ricreare due liste, altrimenti per ogni layer dei plot usa la funzione
plot_dispersion_dc <- wrap_plots(rev(plot_dispersion),ncol = 1)
```
```{r plot dispersion, fig.ext='png', message=FALSE, warning=FALSE, include=FALSE}
dispersione <- filter(rast_df_br,str_detect(data_frame,"dispersione")) #estraggo i dati di dispersione
dispersione$data_frame <- as.factor(str_extract(dispersione$data_frame,"([^\\-]+)(?=\\.\\w+$)")) #rinomino i fattori
ridge_dispersion <- ggplot(dispersione,aes(x=value,y=data_frame,height = stat(count),fill=stat(x)))+ #stat(count) per avere la dimensione legata al numero di celle
geom_density_ridges_gradient(stat = "density",rel_min_height = 0.01,scale=0.9)+ #scale=0.9 per non far toccare le curve
scale_fill_gradientn(colours = brewer.pal(7, "RdYlGn"),
trans = "log10",
label = function(x) sprintf("%.2f", x), #questo per le etichette della legenda
breaks=trans_breaks('log10', function(x) 10^x),
# labels=trans_format('log10', math_format(10^.x)),
limits=c(1, 100) #per avere le legende dei vari plot confrontabili
) +
scale_x_continuous(
trans='log10'
,breaks=trans_breaks('log10', function(x) 10^x),
labels=trans_format('log10', math_format(10^.x))
) +
scale_y_discrete(expand = c(0.01, 0)) +
annotation_logticks(sides = "b")+
labs(x=NULL,y = NULL) +
theme(legend.key.height = unit(50, 'points'),
legend.position = "left")
```
```{r heatmap mortalita, message=FALSE, include=FALSE}
mortalita_df <- rast_df[grepl("mortalita", names(rast_df))] #estraggo solo i dataframe della mortalità
fun_plot_mortality <- function(x) ggplot(x,aes(x=x,y=y))+ #funzione per l'heatmap
geom_tile(aes(fill=value))+
scale_fill_gradientn(colours = rev(brewer.pal(7, "RdYlGn")),
trans='log10',
label = function(x) formatC(x, format = "e", digits = 0),
limits=c(0.0000000000001, 100) #per avere le legende dei vari plot confrontabili
) +
theme(axis.text.x.bottom = element_blank(),
axis.ticks.x=element_blank(),
axis.text.y=element_blank(),
axis.ticks.y=element_blank(),
legend.position = "none",
legend.title = element_blank()
)+
labs(x = NULL, y = NULL) + #tolgo i label agli assi
coord_fixed() #per avere un quadrato sempre
plot_mortality <- lapply(mortalita_df,fun_plot_mortality)
#nomi_plot <- as.list(str_extract(names(mortalita_df),"([^\\-]+)(?=\\.\\w+$)")) #per aggiungere i titoli per ogni plot faccio una lista dei nomi dei raster-dataframe
#plot_mortality <- mapply(plot=plot_mortality,nomi=nomi_plot,fun_nomi_plot,SIMPLIFY = F) #ho due liste, una dei plot e una dei nomi, applico una funzione multipla, simplify=F serve a ricreare due liste, altrimenti per ogni layer dei plot usa la funzione
plot_mortality_dc <- wrap_plots(rev(plot_mortality),ncol = 1)
```
```{r plot mortality, fig.ext='png', message=FALSE, warning=FALSE, include=FALSE}
mortalita <- filter(rast_df_br,str_detect(data_frame,"mortalita")) #estraggo i dati di mortalita
mortalita$data_frame <- as.factor(str_extract(mortalita$data_frame,"([^\\-]+)(?=\\.\\w+$)")) #rinomino i fattori
ridge_mortality <- ggplot(mortalita,aes(x=value,y=data_frame,height = stat(count),fill=stat(x)))+
geom_density_ridges_gradient(stat = "density",rel_min_height = 0.01,scale=0.9)+
scale_fill_gradientn(colours = rev(brewer.pal(7, "RdYlGn")),
trans = "log10",
# label = function(x) sprintf("%.2f", x),
breaks=trans_breaks('log10', function(x) 10^x),
labels=trans_format('log10', math_format(10^.x)),
limits=c(0.0000000000001, 100) #per avere le legende dei vari plot confrontabili
) +
scale_x_continuous(
trans='log10'
,breaks=trans_breaks('log10', function(x) 10^x),
labels=trans_format('log10', math_format(10^.x))
) +
scale_y_discrete(expand = c(0.01, 0)) +
annotation_logticks(sides = "b")+
labs(x=NULL,y = NULL) +
theme(legend.key.height = unit(50, 'points'),
legend.position = "left")
```
#### Foto aeree 2020 AGREA NIR
Queste sono le foto aeree usate per la fotointerpretazione, tratte da [geoportale Emilia-Romagna \> wms \> ortofoto \> AGREA 2020 NIR](https://geoportale.regione.emilia-romagna.it/servizi/servizi-ogc/elenco-capabilities-dei-servizi-wms/cartografia-di-base/service-35). La fotointepretazione è stata supportata da:
- elaborazioni LIDAR risoluzione 1 metro del 2008-2013 quando presenti per l'area di studio, per ricavare la posizione delle arboree e le loro altezze.
- dati AGREA 2021 e 2022, disponibili qui [agreagestione file > AppezzAziendaGrafici](https://agreagestione.regione.emilia-romagna.it/agrea-file/AppezzAziendaGrafici/)
```{r plot foto aeree, echo=FALSE, fig.ext='png',fig.height=6 ,fig.width=16, message=FALSE, warning=FALSE}
do.call(grid.arrange,c(plot_foto_aeree, ncol = (length(raster_files)/(length(raster_dir)-1)),top = "Foto aeree 2020 NIR"))
```
#### Metriche di gradiente di paesaggio
Si tratta di metriche dell'intero paesaggio, senza riferimento alla zona di lancio o presenza del parassitoide.
##### Resistenza
Metrica di quanto è resistente al movimento il paesaggio.
- 1= resistenza totale
- 0= resistenza nulla
```{r plot heatmap resistenza, echo=FALSE, fig.ext='png', fig.width=8 ,fig.height=15, message=FALSE, warning=FALSE}
plot_resistance_dc | ridge_resistance + plot_annotation(title="resistenza")
```
##### Assorbanza
Assorbimento significa impossibilità di movimento, o incapacità di tornare indietro. I valori sono vicini a 1 o vicini a 0.
- 1= assorbanza totale. Nella carta le aree con assorbanza totale sono visualizzate come trasparenti.
- 0= assorbanza nulla, ovvero nessun ulteriore fattore di impossibilità rispetto alla resistenza
```{r plot heatmap assorbimento, echo=FALSE, fig.ext='png', fig.width=8 ,fig.height=15, message=FALSE, warning=FALSE}
plot_absorbance_dc | ridge_assorbance + plot_annotation(title="assorbimento")
```
#### Metriche SAMC (Spatial Absorbing Markov Chain) @samc
Sono le metriche calcolate dal punto di lancio della vespa samurai, che è il centro del raster.
##### Dispersione
Con dispersione si intende la probabilità che la cella sarà visitata almeno una volta dalle 90 femmine lanciate. Non si prendono in considerazione le generazioni successive e il fattore vento.
Il grafico indica quindi la capacità di dispersione delle 90 femmine in quel paesaggio per quel punto di lancio. Nelle mappe sono visualizzati con un gradiente i dati da 1 a 100 per cento. Le zone grigie hanno una probabilità di dispersione inferiore a 1 per cento.
```{r plot heatmap dispersione, echo=FALSE, fig.ext='png',fig.width=8 ,fig.height=15, message=FALSE, warning=FALSE}
plot_dispersion_dc | ridge_dispersion + plot_annotation(title="dispersione")
```
##### Mortalità
Con Mortalità si intende la probabilità di mortalità, che è data dall'insieme dei fattori di resistenza e assorbimento nella cella visitata per quel punto di lancio, secondo la dispersione calcolata. La mortalità è maggiore in vicinanza delle aree di assorbimento, che però per il modello non vengono raggiunte e superate, quindi al loro interno non c'è mortalità. Maggiore è la probabilità di dispersione in una cella maggiore è la probabilità di mortalità.
Nella mappa con scala logaritmica si visualizza la mappa anche nelle zone di dispersione trascurabile (<1%) per agevolare la comprensione del concetto di mortalità
```{r plot heatmap mortalita, echo=FALSE, fig.ext='png',fig.width=8 ,fig.height=15, message=FALSE, warning=FALSE}
plot_mortality_dc | ridge_mortality + plot_annotation(title="mortalità")
```