Hi, First of all I would like to congratulate you for the development of mapview, that now I always use for visualizing spatial data in r, also thanks to the training that I received at the last summer school at the wagening.

Since I upgraded to r 4.0.5, I can no longer render vector data (shapefile) in mackdown pdf document using the mapview function. Instead, I can render raster data in the same document.

I will attach a pdf document with the content of the problem (pp 2-4).

May 04 '21 19:05 agronomofiorentini

Hello, I would like to echo @agronomofiorentini in thanking you for the development of mapview. It is part of my weekly workflows, and I rely on its ability to quickly produce beautiful, portable data visualizations.

I am having the same issue as described by @agronomofiorentini. I am using the mapshot() function to create .html maps from sf objects. Since upgrading to R 4.0.5, the maps produced do not have geometries. The geometries are present in the mapview object when visualized in the Rstudio viewer, but when the mapview object is an argument to mapshot() the resulting html file has the correct extent, correct legend, but no vector features.

May 04 '21 21:05 chrismallon

Can both of you please provide your sessionInfor() after loading mapview?

May 05 '21 06:05 tim-salabim

This is my sessioninfo.

`R version 4.0.5 (2021-03-31) Platform: x86_64-w64-mingw32/x64 (64-bit) Running under: Windows 10 x64 (build 19042)

Matrix products: default

locale: [1] LC_COLLATE=Italian_Italy.1252 LC_CTYPE=Italian_Italy.1252 LC_MONETARY=Italian_Italy.1252 LC_NUMERIC=C LC_TIME=Italian_Italy.1252

attached base packages: [1] stats graphics grDevices utils datasets methods base

other attached packages: [1] mapview_2.9.0

loaded via a namespace (and not attached): [1] Rcpp_1.0.6 pillar_1.6.0 compiler_4.0.5 class_7.3-18 base64enc_0.1-3 tools_4.0.5 digest_0.6.27 satellite_1.0.2
[9] evaluate_0.14 lifecycle_1.0.0 tibble_3.1.1 lattice_0.20-41 pkgconfig_2.0.3 png_0.1-7 rlang_0.4.11 DBI_1.1.1
[17] crosstalk_1.1.1 yaml_2.2.1 xfun_0.22 e1071_1.7-6 dplyr_1.0.5 knitr_1.33 raster_3.4-5 generics_0.1.0
[25] vctrs_0.3.7 htmlwidgets_1.5.3 webshot_0.5.2 classInt_0.4-3 stats4_4.0.5 leaflet_2.0.4.1 DT_0.18 grid_4.0.5
[33] tidyselect_1.1.1 glue_1.4.2 sf_0.9-8 R6_2.5.0 fansi_0.4.2 rmarkdown_2.7 sp_1.4-5 purrr_0.3.4
[41] magrittr_2.0.1 units_0.7-1 scales_1.1.1 codetools_0.2-18 ellipsis_0.3.1 htmltools_0.5.1.1 assertthat_0.2.1 colorspace_2.0-0
[49] KernSmooth_2.23-18 utf8_1.2.1 proxy_0.4-25 munsell_0.5.0 leafem_0.1.3 crayon_1.4.1`

May 05 '21 11:05 agronomofiorentini

Here you are:

R version 4.0.5 (2021-03-31) Platform: x86_64-w64-mingw32/x64 (64-bit) Running under: Windows 10 x64 (build 18363)

Matrix products: default

locale: [1] LC_COLLATE=English_United States.1252 [2] LC_CTYPE=English_United States.1252
[3] LC_MONETARY=English_United States.1252 [4] LC_NUMERIC=C
[5] LC_TIME=English_United States.1252

attached base packages: [1] stats graphics grDevices [4] utils datasets methods
[7] base

other attached packages: [1] reproducible_1.1.1.9007 [2] rmapshaper_0.4.4
[3] mapview_2.9.0
[4] htmltools_0.5.0
[5] viridis_0.5.1
[6] viridisLite_0.3.0
[7] leaflet_2.0.4.1
[8] rlang_0.4.10
[9] tidyr_1.1.3
[10] sp_1.4-2
[11] sf_0.9-8
[12] dplyr_1.0.5

loaded via a namespace (and not attached): [1] httr_1.4.2 bit64_0.9-7 jsonlite_1.7.2 geojsonlint_0.4.0 assertthat_0.2.1 stats4_4.0.5 fpCompare_0.2.3
[8] blob_1.2.1 yaml_2.2.1 remotes_2.1.1 backports_1.1.7
[12] pillar_1.6.0
[13] RSQLite_2.2.0
[14] lattice_0.20-41
[15] glue_1.4.1
[16] quickPlot_0.1.6
[17] digest_0.6.25
[18] RColorBrewer_1.1-2
[19] stringfish_0.12.1
[20] colorspace_1.4-1
[21] leaflet.providers_1.9.0 [22] pkgconfig_2.0.3
[23] httpcode_0.3.0
[24] raster_3.3-6
[25] purrr_0.3.4
[26] scales_1.1.1
[27] webshot_0.5.2
[28] satellite_1.0.2
[29] RApiSerialize_0.1.0
[30] tibble_3.1.1
[31] farver_2.0.3
[32] generics_0.0.2
[33] ggplot2_3.3.3
[34] ellipsis_0.3.1
[35] geojsonio_0.9.4
[36] lazyeval_0.2.2
[37] cli_2.5.0
[38] magrittr_2.0.1
[39] crayon_1.4.1
[40] memoise_1.1.0
[41] maptools_1.0-1
[42] fansi_0.4.1
[43] foreign_0.8-81
[44] class_7.3-18
[45] geojsonsf_2.0.1
[46] tools_4.0.5
[47] data.table_1.12.8
[48] gridBase_0.4-7
[49] lifecycle_1.0.0
[50] geojson_0.3.4
[51] V8_3.4.0
[52] munsell_0.5.0
[53] Require_0.0.4.9002
[54] compiler_4.0.5
[55] e1071_1.7-3
[56] qs_0.22.1
[57] jqr_1.2.0
[58] classInt_0.4-3
[59] units_0.6-7
[60] grid_4.0.5
[61] rstudioapi_0.13
[62] htmlwidgets_1.5.1
[63] igraph_1.2.5
[64] crosstalk_1.1.0.1
[65] leafem_0.1.3
[66] base64enc_0.1-3
[67] gtable_0.3.0
[68] codetools_0.2-18
[69] DBI_1.1.0
[70] jsonvalidate_1.1.0
[71] curl_4.3
[72] R6_2.4.1
[73] gridExtra_2.3
[74] rgdal_1.5-12
[75] rgeos_0.5-3
[76] bit_1.1-15.2
[77] utf8_1.1.4
[78] KernSmooth_2.23-18
[79] crul_1.1.0
[80] Rcpp_1.0.4.6
[81] vctrs_0.3.8
[82] png_0.1-7
[83] tidyselect_1.1.0

May 06 '21 21:05 chrismallon

Can you please remotes::install_github("r-spatial/mapview") and try if that solves it?

May 07 '21 06:05 tim-salabim

Yes, I have reinstalled it but the problem is not solved.

PES_GIOVANNI_3629.pdf

I'm sorry, I hope this doesn't take up a lot of your time. Thanks as always

May 07 '21 06:05 agronomofiorentini

It's fine. Does it work if you explicitly set

mapviewOptions(fgb = FALSE)

before creating the map? I.e. somewhere at the beginning of your script?

May 07 '21 06:05 tim-salabim

Now it is working well.

Again, as always, thanks.

May 07 '21 07:05 agronomofiorentini

@agronomofiorentini do you maybe have a reproducible example of creating a pdf that fails to show the vector geometries without setting mapviewOptions(fgb = FALSE) so that I can investigate how to avoid having to set it explicitly when rendering to pdf?

May 07 '21 09:05 tim-salabim

well it is difficult to share it here, if you want i can share a folder in google drive with all the materials to reproduce exactly what i showed you.

Because to reproduce it you need several files.

May 07 '21 09:05 agronomofiorentini

I was more thinking of a snippet that uses built-in data (e.g. franconia)

May 07 '21 09:05 tim-salabim

I don't know how to do it, sorry.

May 07 '21 09:05 agronomofiorentini

Ok, can you at least post the code? I can easily substitute the relevant bits myself. I just need to know how you specify the markdown

May 07 '21 10:05 tim-salabim

ok

title: "Consorzio Di Bonifica del Nord Sardegna" subtitle: "Report di analisi dell'NDVI" author: "Ernesto Usai (1) & Luigi Ledda (2) & Marco Fiorentini (3)" date: "Maggio 2021" output: pdf_document params: sito: sito DTM: DTM Prima: Prima Prossimo: Prossimo Dopo_30gg: Dopo_30gg dsn: dsn layer: layer dati_danno: dati_danno date: date

knitr::opts_chunk$set(echo = FALSE,
                      cache = FALSE,
                      message = FALSE,
                      warning = FALSE)

library(rgdal)
library(raster)
library(mapview)
library(gt)
library(ggplot2)
library(reshape)
library(emmeans)
library(multcomp)
library(NbClust)
library(rgeos)
library(tidyverse)
library(sp)
library(ggspatial)
library(ggpubr)
library(viridis)

mapviewOptions(fgb=FALSE)
RdYlGr = colorRampPalette(c('Red', 'Yellow', "Green"))

1 - Dottore Agronomo, libero professionista.

2 - Professore Ordinario di Produzione erbacce alimentari presso il Dipartimento di Agraria dell'Università Politecnica delle Marche.

3 - Dottore Agronomo, Ricercatore presso il Dipartimento di Agraria dell'Università Politecnica delle Marche.

$$\[0.1in]$$

Dati richiedente

names(dati_danno)<-c("Richiedente", "Protocollo", "Data evento", "Data protocollo", "Danno")
gt(head(dati_danno,1))

$$\[0.1in]$$

Il presente report si compone di due parti.

Nella prima parte sono raccolti i dati relativi alla richiesta di risarcimento danni formulata dal richiedente, con ubicazione spaziale delle particelle catastali oggetto di richiesta su base foto aerea e indicazione della morfologia del terreno mediante utilizzo di un Modello Digitale del Terreno (DTM) con passo a 1 m.

Nella seconda parte è possibile visualizzare e analizzare le immagini satellitari messe a disposizione gratuitamente dall' European Space Agency (ESA). Le immagini hanno le seguenti caratteristiche:

Immagini appartenenti al Long-Term-Archive (LTA) del sistema Copernicus.
Immagini acquisite con i satelliti Sentinel 2A & 2B.
Immagini riferenti al codice 32TML.
Le immagini scaricate ed elaborate hanno come copertura nuvolosa al massimo del 30%.
Sistema geografico di riferimento WGS 84 - UTM 32 N.
Risoluzione spaziale di 10 x 10 m.

Le immagini satellitari consentono di calcolare l'indice di vegetazione NDVI, utilizzato allo scopo di verificare se a seguito di un evento (esondazione di canali e/o rottura di condotta) si sia avuta una tangibile variazione dell'indice di vegetazione e, in quel caso, su quale superficie dell'azienda.

Si rimanda alla RELAZIONE AGRONOMICA per l'analisi del metodo utilizzato.

\newpage

Sito in analisi

Il sito in analisi è il seguente.

mapview(sito, map.types = c("Esri.WorldImagery"), popup=FALSE)

Dati catastali del sito

I dati catastali del sito in esamina sono riportati di seguito.

Catasto<-sito$CATASTO
Comune<-sito$Comune
Foglio<-sito$FOGLIO
Numero<-sito$NUMERO
Area<-sito$Area
Qualità<-sito$Qualita
Classe<-sito$CLASSE
Reddito_agrario<-sito$RED_AGRARI
Reddito_domenicale<-sito$RED_DOMENI
Cognome<-sito$COGNOME
Nome<-sito$NOME

dati_catastali_1<-data.frame(Catasto, Comune, Foglio, Numero, Area)

dati_catastali_2<-data.frame(Qualità, Classe, Reddito_agrario, Reddito_domenicale)

names(dati_catastali_1)<-c("Catasto*", "Comune", "Foglio", "Particella", "Area (m2)")
gt(dati_catastali_1)

*T=Terreni o F=Fabbricati.

names(dati_catastali_2)<-c("Qualita'", "Classe", "Reddito agrario", "Reddito domenicale")
gt(dati_catastali_2)

\newpage

Modello Digitale del Terreno

IL Modello Digitale del Terreno è utile per rappresentare la morfologia del sito in esame e l'andamento altimetrico delle quote.

La risoluzione spaziale del Modello Digitale del Terreno è di 1x1 m.

Base DTM fornita dal Consorzio di Bonifica del Nord Sardegna.

clip<-function(raster,shape) {
  a1_crop<-crop(raster,shape)
  step1<-rasterize(shape,a1_crop)
  a1_crop*step1}
DTM<-clip(DTM, sito)
names(DTM)<-c("DTM")
mapview(DTM, map.types = c("Esri.WorldImagery"), trim=F)

\newpage

Download, Calcolo e Visualizzazione Indice NDVI

In questa sezione si riporta l'analisi dell'indice NDVI calcolato per il sito:

Prima dell'evento
Subito dopo l'evento
Dopo circa 30 giorni dalla data dell'evento

Le date delle immagini satellitari utilizzate sono le seguenti.

names(date)<-c("Prima dell'evento", "Data dell'evento", "Subito dopo l'evento", "~Dopo 30 gg dall'evento")
gt(date)

Le immagini hanno la medesima scala di rappresentazione per rendere omogenea la comparazione.

#Prima<-clip(Prima, sito)
mapview(Prima, map.types = c("Esri.WorldImagery"), at = seq(0, 1, 0.1), col.regions=RdYlGr)

#Dopo<-clip(Dopo, sito)
mapview(Prossimo, map.types = c("Esri.WorldImagery"), at = seq(0, 1, 0.1), col.regions=RdYlGr)

#Dopo<-clip(Dopo, sito)
mapview(Dopo_30gg, map.types = c("Esri.WorldImagery"), at = seq(0, 1, 0.1), col.regions=RdYlGr)

\newpage

Istogramma di variazione dell'indice

Di seguito è possibile visualizzare l'istogramma del valore medio dell'indice NDVI:

Prima dell'evento
Subito dopo dell'evento.
Dopo circa 30 giorni dall'evento.

before_data_frame<-as.data.frame(Prima)
before_data_frame<-na.omit(before_data_frame)
before_data_mean<-mean(before_data_frame$layer)

prossimo_data_frame<-as.data.frame(Prossimo)
prossimo_data_frame<-na.omit(prossimo_data_frame)
prossimo_data_mean<-mean(prossimo_data_frame$layer)

dopo_data_frame<-as.data.frame(Dopo_30gg)
dopo_data_frame<-na.omit(dopo_data_frame)
dopo_data_mean<-mean(dopo_data_frame$layer)

data_frame_anova<-data.frame(before_data_frame, prossimo_data_frame, dopo_data_frame)
names(data_frame_anova)<-c("Prima", "Prossimo", "~Dopo 30 gg")
data_frame_anova<-melt(data_frame_anova)

data_frame<-data.frame(before_data_mean, prossimo_data_mean, dopo_data_mean)
names(data_frame)<-c("Prima", "Prossimo", "~Dopo 30 gg")
data_frame<-melt(data_frame)

ggplot(data_frame) +
  aes(x = variable, fill = variable, weight = value) +
  geom_bar() +
  scale_fill_viridis_d(option = "viridis") +
  labs(x = "Data", y = "NDVI", title = "NDVI", fill = "Data") +
  theme_gray()

\newpage

Anova

Per poter stimare se, a seguito dell'evento, si è avuta una differenza significativa nel valore dell’indice NDVI, tale cioè da compromettere lo sviluppo della coltura, si è proceduto ad eseguire un’analisi statistica con metodo ANOVA mettendo a confronto gli indici calcolati prima dell'evento, subito dopo l'evento e dopo 30 giorni dall'evento.

L'analisi statistica che si esegue di seguito è composta dai seguenti step:

Creazione di un modello statistico che considera il tempo come fattore sperimentale.
Valutazione se il modello creato soddisfa le 3 assunzioni principali per l'applicazione dell'ANOVA.
Applicazione dell'ANOVA.
Solo se l'ANOVA manifesta una differenza significativa, andiamo ad indagare l'esistenza di differenze significative fra le date riportate.

onefactormodel<-lm(value ~ variable,
                   data=na.omit(data_frame_anova))

anova(onefactormodel)

emmeans_bonf<-emmeans(onefactormodel, ~variable)
cld(emmeans_bonf, reversed = TRUE, Letters=LETTERS, adjust="bonf")

\newpage

Analisi cluster

In questa sezione si riportano i risultati di un'analisi cluster che consente di evidenziare se a seguito dell'evento, si è generata la presenza di zone non appartenenti al medesimo gruppo.

All'interno dei gruppi i dati non hanno alcuna differenza statistica, invece fra diversi gruppi i dati hanno una differenza statistica.

unique_stack<-raster::stack(Prima, Prossimo, Dopo_30gg)
data_frame_monte_xy <- rasterToPoints(unique_stack, spatial=TRUE)
data_frame_monte <- as.data.frame(data_frame_monte_xy)
data_frame_monte<-data_frame_monte[,1:3]
data_frame_monte<-data.frame(data_frame_monte)

cluster<-NbClust(data = data_frame_monte[,-2], 
                 distance = "euclidean",
                 min.nc = 2, 
                 max.nc = 10,
                 index = "gap",
                 method = "kmeans")

Di seguito viene visualizzata l'area per ogni singola zona evidenziata dall'analisi precedente.

xy <- as.data.frame(coordinates(data_frame_monte_xy))
ZMmap_clust <- data.frame(cluster$Best.partition, xy)
ZMmap_d <- data.frame(cluster$Best.partition, xy)
coordinates(ZMmap_d) <- ~ x + y
crs.geo <- CRS("+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs")
proj4string(ZMmap_d) <- crs.geo
gridded(ZMmap_d) <- TRUE
ZMmap_b <- raster(ZMmap_d)
Zone_a <- rasterToPolygons(ZMmap_b)
Zone<-ZMmap_b

max<-Zone %>% 
  as.data.frame() %>% 
  distinct(cluster.Best.partition) %>% 
  na.omit() %>% 
  max()

min<-Zone %>% 
  as.data.frame() %>% 
  distinct(cluster.Best.partition) %>% 
  na.omit() %>% 
  min()

mapview(ZMmap_b, map.types = c("Esri.WorldImagery"), popup=FALSE, at = seq(0, max, 1), col.regions=viridis)

writeOGR(obj = Zone_a,
         dsn = dsn,
         layer = layer,
         driver = "ESRI Shapefile",
         overwrite_layer = T)

Di seguito è riportato invece i grafici dell'indice NDVI suddiviso per le zone generate dall'analisi cluster Subito prima e dopo del danno e dopo 30 gg dal danno.

I grafici possiedono la medesima scala e ciò consente di eseguire una comparazione del valore dell'indice nelle date qua riportate.

DTM_visualization<-data.frame(cluster$Best.partition, data_frame_monte)
DTM_visualization$cluster.Best.partition<-as.character(DTM_visualization$cluster.Best.partition)

Prima<-ggplot(DTM_visualization) +
  aes(x = reorder(cluster.Best.partition, layer.1, na.rm=T), y = layer.1, fill = cluster.Best.partition) +
  geom_boxplot() +
  scale_fill_viridis_d(option = "viridis") +
  labs(x = "Zone", y = "NDVI", title = "Prima", subtitle = "Prima del danno", fill = "Zone") +
  ylim(0L, 1L) +
  theme_gray()

Prossimo<-ggplot(DTM_visualization) +
  aes(x = reorder(cluster.Best.partition, layer.2, na.rm=T), y = layer.2, fill = cluster.Best.partition) +
  geom_boxplot() +
  scale_fill_viridis_d(option = "viridis") +
  labs(x = "Zone", y = "NDVI", title = "Dopo", subtitle = "Subito dopo del danno", fill = "Zone") +
  ylim(0L, 1L) +
  theme_gray()

Dopo<-ggplot(DTM_visualization) +
  aes(x = reorder(cluster.Best.partition, layer.2, na.rm=T), y = layer.3, fill = cluster.Best.partition) +
  geom_boxplot() +
  scale_fill_viridis_d(option = "viridis") +
  labs(x = "Zone", y = "NDVI", title = "Dopo ~30gg", subtitle = "Dopo ~30gg dal danno", fill = "Zone") +
  ylim(0L, 1L) +
  theme_gray()

ggarrange(Prima, Prossimo, Dopo, ncol=3, common.legend = TRUE)

Di seguito invece viene riportata l'area per ogni singola zona evidenziata dall'analisi precedente.

n_clust<-unique(ZMmap_clust$cluster.Best.partition)
n_clust<-as.data.frame(n_clust)
n_clust<-nrow(n_clust)

if (n_clust== 2) {
  
  Zone_1 <- subset(Zone_a, cluster.Best.partition == "1")
  Zone_2 <- subset(Zone_a, cluster.Best.partition == "2")
  Zone_1<-gArea(Zone_1)
  Zone_2<-gArea(Zone_2)
  zone<-data.frame(Zone_1, Zone_2)
  names(zone)<-c("Area Zona 1 (m2)", "Area Zona 2 (m2)")
  gt(zone)
  
} else if (n_clust== 3) {
  
  Zone_1 <- subset(Zone_a, cluster.Best.partition == "1")
  Zone_2 <- subset(Zone_a, cluster.Best.partition == "2")
  Zone_3 <- subset(Zone_a, cluster.Best.partition == "3")
  Zone_1<-gArea(Zone_1)
  Zone_2<-gArea(Zone_2)
  Zone_3<-gArea(Zone_3)
  zone<-data.frame(Zone_1, Zone_2, Zone_3)
  names(zone)<-c("Area Zona 1 (m2)", "Area Zona 2 (m2)", "Area Zona 3 (m2)")
  gt(zone)
  
} else if (n_clust== 4) {
  
  Zone_1 <- subset(Zone_a, cluster.Best.partition == "1")
  Zone_2 <- subset(Zone_a, cluster.Best.partition == "2")
  Zone_3 <- subset(Zone_a, cluster.Best.partition == "3")
  Zone_4 <- subset(Zone_a, cluster.Best.partition == "4")
  Zone_1<-gArea(Zone_1)
  Zone_2<-gArea(Zone_2)
  Zone_3<-gArea(Zone_3)
  Zone_4<-gArea(Zone_4)
  zone<-data.frame(Zone_1, Zone_2, Zone_3, Zone_4)
  names(zone)<-c("Area Zona 1 (m2)", "Area Zona 2 (m2)", "Area Zona 3 (m2)", "Area Zona 4 (m2)")
  gt(zone)
  
} else if (n_clust== 5) {
  
  Zone_1 <- subset(Zone_a, cluster.Best.partition == "1")
  Zone_2 <- subset(Zone_a, cluster.Best.partition == "2")
  Zone_3 <- subset(Zone_a, cluster.Best.partition == "3")
  Zone_4 <- subset(Zone_a, cluster.Best.partition == "4")
  Zone_5 <- subset(Zone_a, cluster.Best.partition == "5")
  Zone_1<-gArea(Zone_1)
  Zone_2<-gArea(Zone_2)
  Zone_3<-gArea(Zone_3)
  Zone_4<-gArea(Zone_4)
  Zone_5<-gArea(Zone_5)
  zone<-data.frame(Zone_1, Zone_2, Zone_3, Zone_4, Zone_5)
  names(zone)<-c("Area Zona 1 (m2)", "Area Zona 2 (m2)", "Area Zona 3 (m2)", "Area Zona 4 (m2)", "Area Zona 5 (m2)")
  gt(zone)
  
} else if (n_clust== 6) {
  
  Zone_1 <- subset(Zone_a, cluster.Best.partition == "1")
  Zone_2 <- subset(Zone_a, cluster.Best.partition == "2")
  Zone_3 <- subset(Zone_a, cluster.Best.partition == "3")
  Zone_4 <- subset(Zone_a, cluster.Best.partition == "4")
  Zone_5 <- subset(Zone_a, cluster.Best.partition == "5")
  Zone_6 <- subset(Zone_a, cluster.Best.partition == "6")
  Zone_1<-gArea(Zone_1)
  Zone_2<-gArea(Zone_2)
  Zone_3<-gArea(Zone_3)
  Zone_4<-gArea(Zone_4)
  Zone_5<-gArea(Zone_5)
  Zone_6<-gArea(Zone_6)
  zone<-data.frame(Zone_1, Zone_2, Zone_3, Zone_4, Zone_5, Zone_6)
  names(zone)<-c("Area Zona 1 (m2)", "Area Zona 2 (m2)", "Area Zona 3 (m2)", "Area Zona 4 (m2)", "Area Zona 5 (m2)", "Area Zona 6 (m2)")
  gt(zone)
  
} else if (n_clust== 7) {
  
  Zone_1 <- subset(Zone_a, cluster.Best.partition == "1")
  Zone_2 <- subset(Zone_a, cluster.Best.partition == "2")
  Zone_3 <- subset(Zone_a, cluster.Best.partition == "3")
  Zone_4 <- subset(Zone_a, cluster.Best.partition == "4")
  Zone_5 <- subset(Zone_a, cluster.Best.partition == "5")
  Zone_6 <- subset(Zone_a, cluster.Best.partition == "6")
  Zone_7 <- subset(Zone_a, cluster.Best.partition == "7")
  Zone_1<-gArea(Zone_1)
  Zone_2<-gArea(Zone_2)
  Zone_3<-gArea(Zone_3)
  Zone_4<-gArea(Zone_4)
  Zone_5<-gArea(Zone_5)
  Zone_6<-gArea(Zone_6)
  Zone_7<-gArea(Zone_7)
  zone<-data.frame(Zone_1, Zone_2, Zone_3, Zone_4, Zone_5, Zone_6, Zone_7)
  names(zone)<-c("Area Zona 1 (m2)", "Area Zona 2 (m2)", "Area Zona 3 (m2)", "Area Zona 4 (m2)", "Area Zona 5 (m2)", "Area Zona 6 (m2)", "Area Zona 7 (m2)")
  gt(zone)
  
} else if (n_clust== 8) {
  
  Zone_1 <- subset(Zone_a, cluster.Best.partition == "1")
  Zone_2 <- subset(Zone_a, cluster.Best.partition == "2")
  Zone_3 <- subset(Zone_a, cluster.Best.partition == "3")
  Zone_4 <- subset(Zone_a, cluster.Best.partition == "4")
  Zone_5 <- subset(Zone_a, cluster.Best.partition == "5")
  Zone_6 <- subset(Zone_a, cluster.Best.partition == "6")
  Zone_7 <- subset(Zone_a, cluster.Best.partition == "7")
  Zone_8 <- subset(Zone_a, cluster.Best.partition == "8")
  Zone_1<-gArea(Zone_1)
  Zone_2<-gArea(Zone_2)
  Zone_3<-gArea(Zone_3)
  Zone_4<-gArea(Zone_4)
  Zone_5<-gArea(Zone_5)
  Zone_6<-gArea(Zone_6)
  Zone_7<-gArea(Zone_7)
  Zone_8<-gArea(Zone_8)
  zone<-data.frame(Zone_1, Zone_2, Zone_3, Zone_4, Zone_5, Zone_6, Zone_7, Zone_8)
  names(zone)<-c("Area Zona 1 (m2)", "Area Zona 2 (m2)", "Area Zona 3 (m2)", "Area Zona 4 (m2)", "Area Zona 5 (m2)", "Area Zona 6 (m2)", "Area Zona 7 (m2)", "Area Zona 8 (m2)")
  gt(zone)
  
} else if (n_clust== 9) {
  
  Zone_1 <- subset(Zone_a, cluster.Best.partition == "1")
  Zone_2 <- subset(Zone_a, cluster.Best.partition == "2")
  Zone_3 <- subset(Zone_a, cluster.Best.partition == "3")
  Zone_4 <- subset(Zone_a, cluster.Best.partition == "4")
  Zone_5 <- subset(Zone_a, cluster.Best.partition == "5")
  Zone_6 <- subset(Zone_a, cluster.Best.partition == "6")
  Zone_7 <- subset(Zone_a, cluster.Best.partition == "7")
  Zone_8 <- subset(Zone_a, cluster.Best.partition == "8")
  Zone_9 <- subset(Zone_a, cluster.Best.partition == "9")
  Zone_1<-gArea(Zone_1)
  Zone_2<-gArea(Zone_2)
  Zone_3<-gArea(Zone_3)
  Zone_4<-gArea(Zone_4)
  Zone_5<-gArea(Zone_5)
  Zone_6<-gArea(Zone_6)
  Zone_7<-gArea(Zone_7)
  Zone_8<-gArea(Zone_8)
  Zone_9<-gArea(Zone_9)
  zone<-data.frame(Zone_1, Zone_2, Zone_3, Zone_4, Zone_5, Zone_6, Zone_7, Zone_8, Zone_9)
  names(zone)<-c("Area Zona 1 (m2)", "Area Zona 2 (m2)", "Area Zona 3 (m2)", "Area Zona 4 (m2)", "Area Zona 5 (m2)", "Area Zona 6 (m2)", "Area Zona 7 (m2)", "Area Zona 8 (m2)", "Area Zona 9 (m2)")
  gt(zone)
  
} else {
  
  Zone_1 <- subset(Zone_a, cluster.Best.partition == "1")
  Zone_2 <- subset(Zone_a, cluster.Best.partition == "2")
  Zone_3 <- subset(Zone_a, cluster.Best.partition == "3")
  Zone_4 <- subset(Zone_a, cluster.Best.partition == "4")
  Zone_5 <- subset(Zone_a, cluster.Best.partition == "5")
  Zone_6 <- subset(Zone_a, cluster.Best.partition == "6")
  Zone_7 <- subset(Zone_a, cluster.Best.partition == "7")
  Zone_8 <- subset(Zone_a, cluster.Best.partition == "8")
  Zone_9 <- subset(Zone_a, cluster.Best.partition == "9")
  Zone_10 <- subset(Zone_a, cluster.Best.partition == "10")
  Zone_1<-gArea(Zone_1)
  Zone_2<-gArea(Zone_2)
  Zone_3<-gArea(Zone_3)
  Zone_4<-gArea(Zone_4)
  Zone_5<-gArea(Zone_5)
  Zone_6<-gArea(Zone_6)
  Zone_7<-gArea(Zone_7)
  Zone_8<-gArea(Zone_8)
  Zone_9<-gArea(Zone_9)
  Zone_10<-gArea(Zone_10)
  zone<-data.frame(Zone_1, Zone_2, Zone_3, Zone_4, Zone_5, Zone_6, Zone_7, Zone_8, Zone_9, Zone_10)
  names(zone)<-c("Area Zona 1 (m2)", "Area Zona 2 (m2)", "Area Zona 3 (m2)", "Area Zona 4 (m2)", "Area Zona 5 (m2)", "Area Zona 6 (m2)", "Area Zona 7 (m2)", "Area Zona 8 (m2)", "Area Zona 9 (m2)", "Area Zona 10 (m2)")
  gt(zone)
  
}

Osservazioni (Eventuale)

\newpage

Differenza statistica fra le zone prima del danno

onefactormodel<-lm(layer.1 ~ cluster.Best.partition,
                   data=na.omit(DTM_visualization))

anova(onefactormodel)

emmeans_bonf<-emmeans(onefactormodel, ~cluster.Best.partition)
cld(emmeans_bonf, reversed = TRUE, Letters=LETTERS, adjust="bonf")

Differenza statistica fra le zone subito dopo del danno

onefactormodel<-lm(layer.2 ~ cluster.Best.partition,
                   data=na.omit(DTM_visualization))

anova(onefactormodel)

emmeans_bonf<-emmeans(onefactormodel, ~cluster.Best.partition)
cld(emmeans_bonf, reversed = TRUE, Letters=LETTERS, adjust="bonf")

\newpage

Differenza statistica fra le zone subito dopo ~30gg dal danno

onefactormodel<-lm(layer.3 ~ cluster.Best.partition,
                   data=na.omit(DTM_visualization))

anova(onefactormodel)

emmeans_bonf<-emmeans(onefactormodel, ~cluster.Best.partition)
cld(emmeans_bonf, reversed = TRUE, Letters=LETTERS, adjust="bonf")

Bibliografia

May 07 '21 10:05 agronomofiorentini

Thanks!

May 07 '21 11:05 tim-salabim

For future reference, here's a minimal reproducible example (note that you need to delete the 4 spaces at the beginning of the lines to get this rendering correctly from RStudio):

    ---
    title: "Some Title"
    subtitle: "A subtitle"
    author: "Max Mustermann"
    date: "Today"
    output: pdf_document
    ---
    
    ```{r setup, include=FALSE}
    knitr::opts_chunk$set(echo = FALSE,
                          cache = FALSE,
                          message = FALSE,
                          warning = FALSE)
    ```
    
    ## Query fgb option and mapview version
    
    ```{r mv, echo=TRUE}
    library(mapview)
    
    mapviewGetOption("fgb")
    packageVersion("mapview")
    
    mapview(franconia, map.types = c("Esri.WorldImagery"), popup=FALSE)
    ```

which renders just fine for me without explicitly setting mapviewOptions(fgb = FALSE). That's the expected behaviour as mapview automatically sets it when in a markdown rendering context.

@agronomofiorentini would you please try this minimal reprex and report back if that does not render properly for you? Might be an issue of operating system (I am on Linux).

May 08 '21 08:05 tim-salabim

Seems to work on MacOS:

May 08 '21 11:05 markolipka

Nice! Thanks!

May 08 '21 11:05 tim-salabim

And on a Windows 10 machine also:

The only difference is the strange grey borders in the webshot maps.

May 08 '21 12:05 markolipka

Sorry tim for the delay of my reply

ciao.pdf

The mapview it still shows the gray board.

But it is working perfectly

Oct 08 '21 15:10 agronomofiorentini

Mapview render raster but not vector file pdf_document rmackdown

Sito in analisi

Dati catastali del sito

Modello Digitale del Terreno

Download, Calcolo e Visualizzazione Indice NDVI

Istogramma di variazione dell'indice

Anova

Analisi cluster

Differenza statistica fra le zone prima del danno

Differenza statistica fra le zone subito dopo del danno

Differenza statistica fra le zone subito dopo ~30gg dal danno

Bibliografia