En este ejercicio, vamos a sacarle partido a los datos que nos hemos molestado en preparar para analizar algunos parámetros de calidad del agua de los ríos del País Vasco para avanzar y profundizar en el manejo y la manipulación de nuestros datos.
Como siempre, preparamos nuestros paquetes a cargar. Por el momento, solo el tidyverse
:
library(tidyverse)
Y a continuación los dos archivos que tenéis disponibles en el ejercicio. Uno, los datos generados en mi bucle. El segundo, un archivo que tenéis disponible también en (http://datos.gob) junto a los datos que nos hemos descargado sobre calidad de agua, que contiene una relación de los nombres de los ríos a los que hacemos referencia y las coordenadas donde los técnicos muestrean. Lo podéis buscar vostros mismos del portal (se llama ‘puntos de muestreo’) o cogerlo del ejercicio.
Este archivo riosPV.csv
ontiene todos los datos muestreados y es fruto del ejercicio guiado de programación.
riosPV<-read_delim("C:/Users/Alejandro/Documents/ISM/riosPV.csv",delim=",")
info_rios<-read_delim("C:/Users/Alejandro/Documents/ISM/Datasets/info_riosPV.csv",delim=";")
head(riosPV)
## # A tibble: 6 x 198
## X1 `Sample Point C~ Year `Estado condici~ `Cianuros total~ Cloruros Cobre
## <dbl> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 AGU126 2015 0.710 5 13.5 1.29
## 2 2 ARA150 2015 0.688 5 25.0 3.92
## 3 3 ART168 2015 0.730 5 33.1 1.63
## 4 4 ART202 2015 0.383 5 108. 5.22
## 5 5 ASU045 2015 0.710 5 20.8 1.12
## 6 6 ASU160 2015 0.681 5 21.4 3.15
## # ... with 191 more variables: `Cobre total` <dbl>, `Coliformes fecales` <dbl>,
## # `Coliformes totales` <dbl>, Conductividad <dbl>, Cromo <dbl>, `Cromo
## # VI+` <dbl>, `DBO5. Demanda Biol<f3>gica de Oxigeno (5 d<ed>as)` <dbl>,
## # `DQO. Demanda Qu<ed>mica de Ox<ed>geno` <dbl>, Dureza <dbl>, `Estreptococos
## # fecales` <dbl>, Fenoles <dbl>, Fluoruros <dbl>, `F<f3>sforo total` <dbl>,
## # Hierro <dbl>, Magnesio <dbl>, Manganeso <dbl>, Mercurio <dbl>,
## # `N<ed>quel` <dbl>, Nitrato <dbl>, Nitrito <dbl>, `Nitr<f3>geno
## # total` <dbl>, `Nitr<f3>geno Total Kjhedahl` <dbl>, Ortofosfato <dbl>,
## # `Ox<ed>geno disuelto` <dbl>, pH <dbl>, Plomo <dbl>, Potasio <dbl>,
## # Selenio <dbl>, Sodio <dbl>, `S<f3>lidos en suspensi<f3>n` <dbl>,
## # Sulfatos <dbl>, `Temperatura del agua` <dbl>, Turbidez <dbl>, Zinc <dbl>,
## # `% Saturaci<f3>n de ox<ed>geno` <dbl>, Amonio <dbl>, `Directiva
## # Vida` <lgl>, `Indice de Calidad General ICG` <dbl>, `Indice Prati` <dbl>,
## # `Estado condiciones generales (<ed>ndice IFQ-R). Ecological Quality
## # Ratio` <dbl>, Amoniaco <dbl>, Alcalinidad <dbl>, `Ars<e9>nico` <dbl>,
## # Bicarbonatos <dbl>, Cadmio <dbl>, Calcio <dbl>, Carbonatos <dbl>,
## # `p.p-DDT` <dbl>, `Parafinas cloradas (Cloroalcanos(C10-C13))` <dbl>,
## # `Pentabromodifenileter (PBDE-100) o 2.2'.4.4'.6-pentabromodifenil
## # eter` <dbl>, `Pentabromodifenileter (PBDE-153) o
## # 2.2'.4.4'.5.5'-hexabromodifenil eter` <dbl>, `Pentabromodifenileter
## # (PBDE-154) o 2.2'.4.4'.5.6'-hexabromodifenil eter` <dbl>,
## # `Pentabromodifenileter (PBDE-28) o 2.4.4'-Tribromodifenil <e9>ter` <dbl>,
## # `Pentabromodifenileter (PBDE-47) o 2.2'.4.4'-tetrabromodfenil eter` <dbl>,
## # `Pentabromodifenileter (PBDE-99)` <dbl>, Pentaclorobenceno <dbl>,
## # Pentaclorofenol <dbl>, Chlorfenvinphos <dbl>, Clorobenceno <dbl>,
## # `Clorpirif<f3>s (Clorpirif<f3>s etil)` <dbl>, `delta-HCH` <dbl>,
## # `Di(2-etilhexil)ftalato (DEHP)` <dbl>, `Diclorobenceno ( suma is<f3>meros
## # orto. meta y para)` <dbl>, Diclorometano <dbl>, Dieldrin <dbl>,
## # `Endosulf<e1>n` <dbl>, `Endosulfan I` <dbl>, `Endosulfan II` <dbl>,
## # `Endosulf<e1>n sulfato` <dbl>, Endrin <dbl>, Etilbenceno <dbl>,
## # Fluoranteno <dbl>, `Gamma-HCH (Lindano)` <dbl>, Heptacloro <dbl>,
## # `Heptacloro ep<f3>xido` <dbl>, Hexaclorobenceno <dbl>,
## # Hexaclorobutadieno <dbl>, `Hexaclorociclohexano (sumatorio
## # m<ed>nimo)` <dbl>, `Indeno(1.2.3-cd)pireno` <dbl>, Isodrin <dbl>,
## # `m+p-xileno` <dbl>, `Metilbenceno (Tolueno)` <dbl>, Metolacloro <dbl>,
## # Naftaleno <dbl>, Nonilfenol <dbl>, Octilfenol <dbl>,
## # `1.1.1-Tricloroetano-Metilcloroformo` <dbl>, `1.2.3-Triclorobenceno` <dbl>,
## # `1.2.4-Triclorobenceno` <dbl>, `1.2-Diclorobenceno=o-Diclorobenceno` <dbl>,
## # `1.2-Dicloroetano` <dbl>, `1.2-Dimetilbenceno (o-Xileno)` <dbl>,
## # `1.3.5-Triclorobenceno` <dbl>, `1.3-Diclorobenceno=m-Diclorobenceno` <dbl>,
## # `1.4-Diclorobenceno=p-Diclorobenceno` <dbl>, Alachlor <dbl>, Aldrin <dbl>,
## # `alfa-HCH` <dbl>, Antraceno <dbl>, `p-p' DDD` <dbl>, ...
head(info_rios)
## # A tibble: 6 x 13
## Code `Sample Point` Type Territory Municipality XETRS89 YETRS89 ZETRS89
## <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl>
## 1 A-BJ~ Olaberria (Ja~ RÍOS Gipuzkoa Irun 596641 4796698 0
## 2 AAR0~ Zamarripa (Ar~ RÍOS Bizkaia Loiu 504575 4794715 NA
## 3 AGU1~ Pandos (Aguer~ RÍOS Bizkaia Trucios-Tur~ 479011 4792248 0
## 4 ARA1~ Egino (Arakil~ RÍOS Araba/Ál~ Asparrena 560046 4746591 0
## 5 ART0~ Iruzubieta (A~ RÍOS Bizkaia Ziortza-Bol~ 538399 4788966 0
## 6 ART1~ Aretxinaga (A~ RÍOS Bizkaia Markina-Xem~ 540903 4791898 NA
## # ... with 5 more variables: `Water range` <chr>, `Water range type` <chr>,
## # Basin <chr>, Section <chr>, Subgroup <chr>
Lo primero de todo es sacarle partido a esta nueva tabla de datos porque nos va a decir entre otras cosas el nombre del río, dato fundamental que hemos obviado a esta hora. Como siempre, el paquete dplyr
de tidyverse tiene las utilidades que necesitamos, en este caso gracias a las funciones de la familia join()
. Lleva un tiempo dominar toda la lógica y las utilidades. Pero en nuestro caso tenemos lo siguiente:
que agregar por la derecha a nuestra tabla principal (los datos) la tabla anexa (nombres, coordenadas) etc.
Reconoceremos qué dato de la segunda tabla tiene que vincularse a la primera por el campo común: Code
y Sample.Point.Code.
Como no se llaman igual, tendremos que renombrar una de las dos columnas. También vais a ver que voy a ir renombrando las columnas con nombres que tienen nombres problemáticos, como los compuestos por más de una palabra
Los códigos de los ríos aparecen 5 veces (tenemos registros de ese punto de 5 años). En la tabla auxiliar una vez. Nuestra función tendrá que devolver todas las combinaciones que encuentra que coinciden los códigos.
La función left_join()
nos vale para esta tarea.
dim(riosPV)
## [1] 670 198
dim(info_rios)
## [1] 203 13
riosPV<-rename(riosPV,"Code"="Sample Point Code")
riosPV<-rename(riosPV,"Condgenerales"="Estado condiciones generales (<ed>ndice IFQ-R)")
riosPV<-left_join(riosPV,info_rios,by=c("Code"))
dim(riosPV)
## [1] 670 210
Ya veremos cómo le podemos sacar partido a estos datos. Pero para empezar, puedo hacerme ya muchas preguntas. Por ejemplo, comprobar si algún dato de calidad ha evolucionado en este tiempo o qué datos son realmente extremos para localizar zonas problemáticas. En general vamos a encontrar una cierta estabilidad en las condiciones, salvo por casos concretos que nos vamos a fijar. Por ejemplo, los niveles inusuales de mercurio de 2015.
ggplot(riosPV,aes(x=as.factor(Year),y=Condgenerales))+geom_boxplot()+
theme_bw()+ylab("Condiciones generales")+xlab("Año")+ggtitle("Evolución general")
## Warning: Removed 113 rows containing non-finite values (stat_boxplot).
riosPV<-rename(riosPV,"Estreptococos"="Estreptococos fecales")
ggplot(riosPV,aes(x=as.factor(Year),y=Estreptococos))+geom_boxplot()+
theme_bw()+ylab("Estreptococos")+xlab("Año")+ggtitle("Evaluación estreptococos fecales")
## Warning: Removed 9 rows containing non-finite values (stat_boxplot).
ggplot(riosPV,aes(x=as.factor(Year),y=Mercurio))+geom_boxplot()+
theme_bw()+ylab("Mercurio")+xlab("Año")+ggtitle("Niveles de mercurio")
## Warning: Removed 5 rows containing non-finite values (stat_boxplot).
Podemos ver los datos del mercurio. En 2015 hay una serie de ríos que tienen unos niveles de mercurio anormalmente altos. Vamos a ver cuáles son:
mercurio<-subset(riosPV,Year==2015 & Mercurio>2)
mercurio
## # A tibble: 23 x 210
## X1 Code Year Condgenerales `Cianuros total~ Cloruros Cobre `Cobre total`
## <dbl> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 2 ARA1~ 2015 0.688 5 25.0 3.92 1.30
## 2 4 ART2~ 2015 0.383 5 108. 5.22 1.72
## 3 6 ASU1~ 2015 0.681 5 21.4 3.15 2.75
## 4 12 BAR1~ 2015 0.749 5 14.6 3.87 1.00
## 5 15 BJA0~ 2015 0.658 5 16.5 12.9 6.27
## 6 18 BUT2~ 2015 0.618 5 27.7 1.53 2.18
## 7 25 DEB4~ 2015 0.620 5 29.7 6.18 2.05
## 8 39 GAL0~ 2015 0.687 5 13.1 9.67 1.39
## 9 41 GOB0~ 2015 0.550 5 36.6 34.1 8.05
## 10 47 IBA1~ 2015 0.602 5 21.7 3.57 1.77
## # ... with 13 more rows, and 202 more variables: `Coliformes fecales` <dbl>,
## # `Coliformes totales` <dbl>, Conductividad <dbl>, Cromo <dbl>, `Cromo
## # VI+` <dbl>, `DBO5. Demanda Biol<f3>gica de Oxigeno (5 d<ed>as)` <dbl>,
## # `DQO. Demanda Qu<ed>mica de Ox<ed>geno` <dbl>, Dureza <dbl>,
## # Estreptococos <dbl>, Fenoles <dbl>, Fluoruros <dbl>, `F<f3>sforo
## # total` <dbl>, Hierro <dbl>, Magnesio <dbl>, Manganeso <dbl>,
## # Mercurio <dbl>, `N<ed>quel` <dbl>, Nitrato <dbl>, Nitrito <dbl>,
## # `Nitr<f3>geno total` <dbl>, `Nitr<f3>geno Total Kjhedahl` <dbl>,
## # Ortofosfato <dbl>, `Ox<ed>geno disuelto` <dbl>, pH <dbl>, Plomo <dbl>,
## # Potasio <dbl>, Selenio <dbl>, Sodio <dbl>, `S<f3>lidos en
## # suspensi<f3>n` <dbl>, Sulfatos <dbl>, `Temperatura del agua` <dbl>,
## # Turbidez <dbl>, Zinc <dbl>, `% Saturaci<f3>n de ox<ed>geno` <dbl>,
## # Amonio <dbl>, `Directiva Vida` <lgl>, `Indice de Calidad General
## # ICG` <dbl>, `Indice Prati` <dbl>, `Estado condiciones generales (<ed>ndice
## # IFQ-R). Ecological Quality Ratio` <dbl>, Amoniaco <dbl>, Alcalinidad <dbl>,
## # `Ars<e9>nico` <dbl>, Bicarbonatos <dbl>, Cadmio <dbl>, Calcio <dbl>,
## # Carbonatos <dbl>, `p.p-DDT` <dbl>, `Parafinas cloradas
## # (Cloroalcanos(C10-C13))` <dbl>, `Pentabromodifenileter (PBDE-100) o
## # 2.2'.4.4'.6-pentabromodifenil eter` <dbl>, `Pentabromodifenileter
## # (PBDE-153) o 2.2'.4.4'.5.5'-hexabromodifenil eter` <dbl>,
## # `Pentabromodifenileter (PBDE-154) o 2.2'.4.4'.5.6'-hexabromodifenil
## # eter` <dbl>, `Pentabromodifenileter (PBDE-28) o 2.4.4'-Tribromodifenil
## # <e9>ter` <dbl>, `Pentabromodifenileter (PBDE-47) o
## # 2.2'.4.4'-tetrabromodfenil eter` <dbl>, `Pentabromodifenileter
## # (PBDE-99)` <dbl>, Pentaclorobenceno <dbl>, Pentaclorofenol <dbl>,
## # Chlorfenvinphos <dbl>, Clorobenceno <dbl>, `Clorpirif<f3>s (Clorpirif<f3>s
## # etil)` <dbl>, `delta-HCH` <dbl>, `Di(2-etilhexil)ftalato (DEHP)` <dbl>,
## # `Diclorobenceno ( suma is<f3>meros orto. meta y para)` <dbl>,
## # Diclorometano <dbl>, Dieldrin <dbl>, `Endosulf<e1>n` <dbl>, `Endosulfan
## # I` <dbl>, `Endosulfan II` <dbl>, `Endosulf<e1>n sulfato` <dbl>,
## # Endrin <dbl>, Etilbenceno <dbl>, Fluoranteno <dbl>, `Gamma-HCH
## # (Lindano)` <dbl>, Heptacloro <dbl>, `Heptacloro ep<f3>xido` <dbl>,
## # Hexaclorobenceno <dbl>, Hexaclorobutadieno <dbl>, `Hexaclorociclohexano
## # (sumatorio m<ed>nimo)` <dbl>, `Indeno(1.2.3-cd)pireno` <dbl>,
## # Isodrin <dbl>, `m+p-xileno` <dbl>, `Metilbenceno (Tolueno)` <dbl>,
## # Metolacloro <dbl>, Naftaleno <dbl>, Nonilfenol <dbl>, Octilfenol <dbl>,
## # `1.1.1-Tricloroetano-Metilcloroformo` <dbl>, `1.2.3-Triclorobenceno` <dbl>,
## # `1.2.4-Triclorobenceno` <dbl>, `1.2-Diclorobenceno=o-Diclorobenceno` <dbl>,
## # `1.2-Dicloroetano` <dbl>, `1.2-Dimetilbenceno (o-Xileno)` <dbl>,
## # `1.3.5-Triclorobenceno` <dbl>, `1.3-Diclorobenceno=m-Diclorobenceno` <dbl>,
## # `1.4-Diclorobenceno=p-Diclorobenceno` <dbl>, Alachlor <dbl>, Aldrin <dbl>,
## # `alfa-HCH` <dbl>, Antraceno <dbl>, `p-p' DDD` <dbl>, `p-p' DDE` <dbl>, ...
Tenemos una gran cantidad de datos, de muchos ríos, muchos entornos que puede ser algo difícil de ver en conjunto. Los análisis multivariantes, que analizan efectivamente estos conjuntos, pueden ser de valiosa ayuda. Vamos a hacer, por ejemplo, un sencillo análisis de componentes principales.
Pero aquí tenemos mezclados los conjuntos de datos de varios años. La información es redundante. Como al principio he visto algunos puntos de peor calidad en los primeros gráficos, me quedo solo con los datos de 2015, para ver si puedo caracterizar algunos puntos.
seleccion_rios<-subset(riosPV,Year==2015)
Muchas de las variables no cuentan con datos completos. Vamos a elegir una selección significativamente más limitada con la función select()
. Voy a elegir parámetros que tengan la información completa y sean relevantes, como los datos de salinidad y minerales, presencia de metales pesados, coliformes.. y al principio del todo, voy a guardar los nombres de los ríos así como la provincia.
seleccion_rios<-seleccion_rios%>%
select("Sample Point","Territory","Cianuros totales","Cloruros","Cobre","Coliformes fecales","Conductividad","Hierro","Plomo","Carbonatos","Calcio","Alcalinidad","Bicarbonatos","Cadmio","pH","Ortofosfato","Nitrato","Magnesio","Estreptococos")%>%na.omit()
head(seleccion_rios)
## # A tibble: 6 x 19
## `Sample Point` Territory `Cianuros total~ Cloruros Cobre `Coliformes fec~
## <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Pandos (Aguer~ Bizkaia 5 13.5 1.29 5300
## 2 Egino (Arakil~ Araba/Ál~ 5 25.0 3.92 1158.
## 3 Ribera (Artib~ Bizkaia 5 33.1 1.63 4630
## 4 Gardotza (Art~ Bizkaia 5 108. 5.22 1127683.
## 5 Zamudio (Asua~ Bizkaia 5 20.8 1.12 2232
## 6 Sangroniz (As~ Bizkaia 5 21.4 3.15 5720
## # ... with 13 more variables: Conductividad <dbl>, Hierro <dbl>, Plomo <dbl>,
## # Carbonatos <dbl>, Calcio <dbl>, Alcalinidad <dbl>, Bicarbonatos <dbl>,
## # Cadmio <dbl>, pH <dbl>, Ortofosfato <dbl>, Nitrato <dbl>, Magnesio <dbl>,
## # Estreptococos <dbl>
Cargamos el paquete Factominer para realizar un PCA algo más visual.
library(FactoMineR)
## Warning: package 'FactoMineR' was built under R version 3.6.3
Realizamos un PCA con los datos numéricos obviando la primera columna con los nombres y también visualizamos la agrupación de los datos. No están particularmente separados, parece que hay una masa de ríos bastante homogénea. Sin embargo hay puntos que se alejan de esa masa siguiendo los ejes del PCA, como el etiquetado con el 90 o el 78 en el eje horizontal, o el 4, o el 41, el 130 o el 103 en el vertical.
En el segundo gráfico que nos muestra la agrupación, importancia y dirección de las variables, vemos que precisamente en ese sentido hay dos agrupaciones de variables: las de metales pesados y bacterias peligrosas en el eje 2 (vertical) y las relacionadas con el incremento de sales de materiales básicos en el eje 1 u horizontal. Esto nos indica que hay algunos ríos con características especiales. ¿Cuáles? Una manera es etiquetar los puntos cambiando los nombres de las filas de mis datos de partida, diciéndole que son los que están en la primera columna, a la que le cambio el nombre para evitar espacios:
seleccion_rios<-rename(seleccion_rios, "Nombre"="Sample Point")
library(ggfortify)
## Warning: package 'ggfortify' was built under R version 3.6.3
library(factoextra)
## Warning: package 'factoextra' was built under R version 3.6.3
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
pca_rios<-PCA(seleccion_rios[3:ncol(seleccion_rios)])
fviz_pca_ind(pca_rios,
label = "ind",
habillage = as.factor(seleccion_rios$Territory), # color by groups
palette = c("#00AFBB", "#E7B800","#FC4E07"),
addEllipses = TRUE # Concentration ellipses
)
rios_contaminados<-seleccion_rios[c(4,128,41,103),]
rios_basicos<-seleccion_rios[c(70,98),]
Parece que los ríos alaveses son los que más se alejan de la masa de ríos, son los que tienen características de mayores concentraciones de sales que alcalinizan, por discurrir por zonas calizas.
Vamos a valernos de herramientas GIS para visualizar estos datos problemáticos. Para avanzar con las herramientas disponibles en R, vamos a ver otro de los paquetes fundamentales, rgdal
, que tendréis que instalar y cargar:
library(rgdal)
## rgdal: version: 1.5-12, (SVN revision 1018)
## Geospatial Data Abstraction Library extensions to R successfully loaded
Mediante la función readOGR
vamos a cargar el shapefile de los límites de las provincias del País Vasco que tenéis disponible en el ejercicio. Yo lo he descomprimido y lo he dejado en una carpeta “Provincia” dentro de mi directorio habitual de datasets.
Euskadi<-readOGR(dsn = "C:/Users/Alejandro/Documents/ISM/Datasets/Provincias/CB_CAPV_5000_ETRS89.shp", stringsAsFactors = F)
## Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS =
## dumpSRS, : Discarded datum European_Terrestrial_Reference_System_1989 in CRS
## definition: +proj=utm +zone=30 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m
## +no_defs
## OGR data source with driver: ESRI Shapefile
## Source: "C:UsersAlejandroDocumentsISMDatasetsProvinciasCB_CAPV_5000_ETRS89.shp", layer: "CB_CAPV_5000_ETRS89"
## with 1 features
## It has 1 fields
Y también el shapefile de las cuencias hidrográficas que delimitan los ríos. Se puede descargar aquí comprimido. También lo he descomprimido y guardado en una carpeta, riosshp.
riosPV_shape<-readOGR(dsn = "C:/Users/Alejandro/Documents/ISM/Datasets/riosshp/CT_0202LRiosCAPV_ETRS89.shp",
stringsAsFactors = F)
## Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS =
## dumpSRS, : Discarded datum European_Terrestrial_Reference_System_1989 in CRS
## definition: +proj=utm +zone=30 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m
## +no_defs
## Warning in readOGR(dsn = "C:/Users/Alejandro/Documents/ISM/Datasets/riosshp/
## CT_0202LRiosCAPV_ETRS89.shp", : Deleted feature IDs: 118, 119
## OGR data source with driver: ESRI Shapefile
## Source: "C:UsersAlejandroDocumentsISMDatasetsriosshpCT_0202LRiosCAPV_ETRS89.shp", layer: "CT_0202LRiosCAPV_ETRS89"
## with 13982 features
## It has 39 fields
## Integer64 fields read as strings: OBJECTID
Ya comentábamos que los shapefiles tienen bastantes características así como estructuras complejas de bases de datos, topología, atributos… EN este caso tenemos un archivo ciertamente pesado, porque recoge con gran detalle todas las cuencas de los ríos, grandes y pequeños. Vamos a simplificarlo y nos vamos a quedar solo con los que tenga jerarquía 1, que recoge los ríos de un tamaño razonablemente grande. Al menos, una operación subset nos permite realizarlo de manera sencilla…aunque es posible que tarde un poco.
riosPV_selection<-subset(riosPV_shape,JERARQUIA==1)
Ahora vamos a mapear todos los puntos que en 2015 nos dieron datos conflictivos de mercurio. En verde vamos a poner todos aquellos puntos de muestreo con datos razonables. En rojo, todos los puntos conflictivos. Pues bien, ggplot nos va a permitir realizar el mapa, aceptando los polígonos que hemos cargado.
El contorno de Euskadi es un polígono: lo dibujaremos con geom_polygon()
; Hay que añadir el argumento group=group
para que se dibuje correctamente.
Para dibujar las líneas de los ríos, geom_path()
.
Y para dibujar los datos de los puntos, con un scatterplot típico como otras veces, nos vale. Primero lo hacemos con todos los datos de muestreo en verde, y luego poniendo otra capa encima solo con los conflictivos, en rojo, tenemos lo que queremos.
ggplot(data=Euskadi)+
geom_polygon(aes(x=long,y=lat,group=group),fill="white",col="black")+
geom_path(data=riosPV_selection,aes(long,lat,group=group),col="blue")+xlab("Longitud")+
geom_point(data=riosPV,aes(x=XETRS89,y=YETRS89),col="green")+
geom_point(data=mercurio,aes(x=XETRS89,y=YETRS89),col="red")+
ylab("Latitud")+theme_bw()+ggtitle("Presencia excesiva de mercurio en el País Vasco")
Este mapa tiende a indicarnos una cierta concentración de los puntos más duros cerca de la Ría de Bilbao, uno de los corazones industriales de la potente industria vasca. Pero no es el único, hay cierta tendencia a concentrar los peores puntos en las desembocaduras al mar. Ahora vamos a visualizar, en vez de puntos rojos, etiquetas sobre el mapa que nos den el nombre de las estaciones contaminadas. Si lo hacemos con ggplot solo, tendremos el problema de que se marcan las etiquetas justo encima de los puntos y suele ser engorroso. Es justo mi propósito: enseñaros una de las extensiones de ggplot que más vais a ver en los manuales, que es ggrepel
y que nos permite desplazar las etiquetas. Tendréis que descargar el paquete si no lo hicisteis en Rmarkdown.
install.packages("ggrepel")
Las funciones de esta extensión son casi idénticas pero añadiendo repel()
al final de la función. Así, geom_label_repel()
es lo mismo que geom_label()
pero evitando poner la etiqueta justo encima de las coordenadas que le estamos indicando. Además hay que usar el argumento label=
para decirle de qué columna de nuestros datos queremos que saque el texto de la etiqueta.
library("ggrepel")
## Warning: package 'ggrepel' was built under R version 3.6.3
ggplot(data=Euskadi)+
geom_polygon(aes(x=long,y=lat,group=group),fill="white",col="black")+
geom_path(data=riosPV_selection,aes(long,lat,group=group),col="blue")+xlab("Longitud")+
geom_point(data=riosPV,aes(x=XETRS89,y=YETRS89),col="green")+
geom_point(data=mercurio,aes(x=XETRS89,y=YETRS89),col="red")+
geom_label_repel(data=mercurio,aes(x=XETRS89,y=YETRS89,label=Basin),size=2.8)+
ylab("Latitud")+theme_bw()
## Regions defined for each Polygons
## Warning: Removed 6 rows containing missing values (geom_point).
Para concluir, vamos a realizar una aplicación GIS con una librería de Javascript llamada leaflet
que nos permite cargar información WMS
o Web Map Service, es decir, servidores GIS que proporcionan información ambiental, de la que tenéis en cualquier visor, y que se han revelado como grandes aliados de los gestores GIS.
Como siempre, primero instalamos leaflet, pero también la libreria sp
si no la tenemos instalada, aunque a estas alturas debería estarlo: install.packages("leaflet","sp")
.
library(leaflet)
## Warning: package 'leaflet' was built under R version 3.6.3
library(sp)
Vamos a realizar un mapa web centrado en la peor medición de mercurio de 2015, alojando mis datos en el objeto mercurio_max. Simplemente seleccionamos aquel dato de mercurio máximo, para comprobar que encuentra en el río butrón o butroi
mercurio_max<-subset(mercurio,Mercurio==max(mercurio$Mercurio))
mercurio_max$Basin
## [1] "Butroe"
A continuación, vamos a extraer las coordenadas, que están en UTM y unidades métricas, para cambiar el sistema de proyeccción espacial y expresarlas en grados:
mercurioUTM<-mercurio_max[203:204]%>%as.data.frame()
colnames(mercurioUTM)<-c("x","y")
Ahora comenzamos la transformación
La función coordinates()
nos permite decirle a R que nos reconozca estos dataframe de dos columnas x e y como coordenadas, ya no como unos puntos cualquiera.
La función CRS()
nos permite crear los parámetros del sistema de proyección UTM en el que se encuentra nuestro punto o puntos. En este caso, decirle que nos saque los parámetros del sistema de proyección denominado EPSGE25830. El dato que nos sale se parece a la cadena de texto que creamos en el objecto projection
que en este caso es un sistema de coordenadas geográficas en grados decimales, y no en metros. Lo pongo de manera diferente para que veais cuál es el resultado de CRS()
.
proj4string()
nos permite asignarle a nuestras nuevas coordenadas geográficas la información de la proyección. Antes solo teníamos el número pero nada decía a qué hacía referencia.
Y finalmente spTransform()
nos permite calcular nuestras coordenadas en metros, a grados decimales:
coordinates(mercurioUTM)<-~x+y
projectionUTM<-CRS('+init=epsg:25830')
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum European_Terrestrial_Reference_System_1989 in CRS definition
projection<-c('+proj=longlat +datum=WGS84 +no_defs')
projection<-
proj4string(mercurioUTM)<-projectionUTM
mercurio_geo<-spTransform(mercurioUTM,projection)%>%as.data.frame()
Las coordenadas en grados es lo que nos permite decirle a las siguientes líneas dónde está centrado mi mapa. Pero, ¿qué mapa?:
En el objeto mapabase me voy a guardar la dirección del servidor de cartografía básica del Instituto Geográfico Nacional. Tenéis a disposición varias decesnas de servidores de ortofotos, cartografía temática y similar.
Inicio la función leaflet()
y con mi operador %>%
voy incorporando elementos:
setView()
nos permite decirle dónde quiero centrar mi mapa y con qué grado de zoom.
addWMSTiles()
nos permite añadir un mapa: quiero que saque la información del objeto mapabase, la capa llamada “IGNBaseOrto”, y que saque el mapa en formato PNG.
Y además le voy a añadir con addCircleMarkers()
un marcador enorme justo donde está la estación de muestreo.
mapabase<-"https://www.ign.es/wms-inspire/ign-base"
leaflet() %>%
setView(lng = mercurio_geo$x, lat = mercurio_geo$y, zoom = 15) %>%
addWMSTiles(
mapabase,
layers = "IGNBaseOrto",
options = WMSTileOptions(format = "image/png", transparent = TRUE)
)%>%
addCircleMarkers(mapabase, lng = mercurio_geo$x, lat = mercurio_geo$y, radius = 10)
En vuestro visor de Rstudio se abrirá una nueva ventana con el mapa. Dependiendo de vuestra máquina y su configuración puede que no ocurra. En ese caso, con el botón que sale para disponer vuestro mapa en el navegador externo, podréis visualizarlo desde allí.
Botón
Vamos a repetir el mapa con los ríos que destacan del análisis multivariante y añadiendo otros elementos. Probad, una vez que abráis el mapa en el navegador, a pasar el cursor por encima de los marcadores (label
) o a pinchar encima (popup
). También la función addTiles()
que nos permite añadir un mapa de OpenStreetMaps.
rios_destacados<-subset(riosPV,Year==2015)
rios_destacados<-rios_destacados[c(4,28,41,103,33,70,98),]
rios_destacados$condicion<-c("Contaminado","Contaminado","Contaminado","Contaminado","Alcalino","Alcalino","Alcalino")
#rios_destacados<-rios_destacados[203:204]%>%as.data.frame()
rios_destacados<-rename(rios_destacados,x=XETRS89,y=YETRS89)
coordinates(rios_destacados)<-~x+y
projectionUTM<-CRS('+init=epsg:25830')
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum European_Terrestrial_Reference_System_1989 in CRS definition
projection<-'+proj=longlat +datum=WGS84 +no_defs'
proj4string(rios_destacados)<-projectionUTM
rios_destacados_geo<-spTransform(rios_destacados,projection)%>%as.data.frame()
Euskadi_geo<-spTransform(Euskadi,projection)
class(rios_destacados_geo)
## [1] "data.frame"
leaflet(Euskadi_geo) %>%
#setView(lng = rios_destacados_geo$x, lat = rios_destacados_geo$y, zoom = 10) %>%
addPolygons(color = "#444444")%>%
addTiles()%>%
addMarkers(rios_destacados_geo,lng=~rios_destacados_geo$x,lat=~rios_destacados_geo$y,
label=~rios_destacados_geo$Basin,popup=~rios_destacados_geo$condicion)