Ejercicio guiado 10.4 (opcional). Análisis de calidad de ríos y GIS

Introducción

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.

Unión de los datos

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:

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

Algunas visualizaciones

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>, ...

Análisis multivariante: caracterización del conjunto

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.

R como SIG

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:

Lectura de datos vectoriales

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.

Mapas con GGPLOT

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").

Creando un servidor SIG

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

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?:

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)