Anàlisi, disseny i implementació d’un quadre de comandament (dashboard) interactiu a la web basat en R-Shiny per l’anàlisi i explotació de la Base de Dades del Mapa Hidrogeològic de Catalunya

Màster Interuniversitari en Estadística i Investigació Operativa UPC-UB Títol: Aplicació interactiva per la consulta i anàlisi de dades hidrogeològiques de Catalunya Autor: Ernest Orriols Director: Ignasi Herms i Isaac Subirana Departament: Institut Cartogràfic i Geològic de Catalunya i Departament de Genètica, Microbiologia i Estadística Universitat: Universitat Politècnica de Catalunya Convocatòria: Juny 2023 : Universitat Politècnica de Catalunya Facultat de Matemàtiques i Estadística Tesi de màster Aplicació interactiva per la consulta i anàlisi de dades hidrogeològiques de Catalunya Ernest Orriols Bertran Directors: Ignasi Herms Isaac Subirana Institut Cartogràfic i Geològic de Catalunya Departament de Genètica, Microbiologia i Estadística Agraïments Vull donar les gràcies i el meu agraïment a la Georgina Arnó, Ignasi Herms i l’Isaac Subirana per ajudar-me a fer aquest treball. Resum Aquesta memòria presenta la implementació d’una eina visual i intuïtiva utilitzant diferents paquets de R com shiny,shinydashboard ileaflet. Aplicant models de classificació no supervisats per la classificació de punts fent ús de l’algoritme de mescles gaussianes del paquet mclust. Prèviament, les dades han estat transformades perquè són dades composicionals. L’aplicació es divideix en cinc apartats en el qual primer hi ha el mapa on es fa la selecció de punts, després la base de dades dels punts seleccionats i finalment l’anàlisi univariant i bivariant i multivariant de les dades seleccionades. En la pàgina multivariant hi ha la funcionalitat més important de l’aplicació, ja que tots els resultats de totes les altres pàgines s’actualitzen d’acord amb les classes o grups trobats a partir del model de mescles gaussianes. Per últim, es fa un cas pràctic que és d’especial interès, perquè s’intenta replicar el que es va fer en un estudi previ realitzat per l’Institut Cartogràfic on s’arriben a conclusions similars. Aquesta aplicació està pública al servidor shinyapps.io amb l’enllaç: https://icgcshiny01022022.shinyapps.io/shinyICGC/ Paraules clau: Shiny, Dashboard, ICGC, Dades composicionals, Clustering Abstract This thesis presents the implementation of a visual and intuitive tool using different R packages such as shiny,shinydashboard and leaflet. Applying unsupervised classification models by point classification using the Gaussian mixture algorithm from the mclust package. Previously, the data has been transformed because they are compositional data. The application is divided into five sections in which first there is the map where the selection of points is done, then the database of the selected points and finally the univariate, bivariate and multivariate analysis of the selected data. On the multivariate page there is the most important functionality of the application, as all results on all other pages are updated according to the classes or groups found from the Gaussian mixture model. Finally, a case is made practical that is of particular interest, because it tries to replicate what was done in a previous study carried out by the Cartographic Institute where similar conclusions are reached. This app is public on shinyapps.io server with the link: https://icgcshiny01022022.shinyapps.io/shinyICGC/ Keywords: Shiny, Dashboard, ICGC, Composicional data, Clustering Índex 1 Introducció 1 2 Objectius 2 3 Mètodes 3 3.1 Paquets de R ............................................. 3 3.1.1 Shiny ............................................. 3 3.1.2 Shinydashboard ........................................ 3 3.1.3 Leaflet ............................................. 3 3.1.4 Mclust ............................................. 4 3.1.5 Plotly ............................................. 4 3.2 Aplicació ............................................... 5 3.2.1 Mapa ............................................. 5 3.2.2 Diagrama de Piper ...................................... 7 3.2.3 Base de dades ......................................... 8 3.2.4 Anàlisi univariant ...................................... 9 3.2.5 Anàlisi bivariant ....................................... 9 3.2.6 Anàlisi multivariant ..................................... 10 3.2.7 Connexió entre les pàgines de l’aplicació .......................... 11 3.3 Eines estadístiques .......................................... 12 3.3.1 Dades composicionals .................................... 12 3.3.2 Model de mescles gaussianes ................................ 12 4 Cas pràctic 14 5 Conclusions 18 6 Annex 19 6.1 Definició de les variables ....................................... 19 6.1.1 Fonts ............................................. 19 6.1.2 Pous .............................................. 21 Bibliografia 23 Índex de figures 1 Pàgina principal mapa. ........................................ 5 2 Exemple d’ús de l’aplicació. ..................................... 6 3 Diagrama de Piper interpretació .................................. 7 4 Diagrama de Piper .......................................... 7 5 Pàgina Base de dades. ........................................ 8 6 Pàgina univariant. .......................................... 9 7 Pàgina bivariant. ........................................... 10 8 Pàgina multivariant. ......................................... 10 9 Pàgina principal mapa actualitzat. ................................. 11 10 Pàgina principal mapa cas pràctic. ................................. 14 11 Pàgina univariant cas pràctic. .................................... 14 12 Pàgina multivariant cas pràctic. ................................... 15 13 Pàgina multivariant cas pràctic. ................................... 15 14 Pàgina principal mapa actualitzat cas pràctic. ........................... 16 Índex de taules 1 Llista de paquets de R utilitzats. .................................. 4 1 Introducció Les aigües subterrànies són un recurs essencial i una peça clau en el cicle de l’aigua i poden ser la solució a les crisis hídriques si es gestionen correctament. Actualment, aporten aproximadament la meitat del volum d’aigua extreta per ús domèstic per la població mundial i al voltant del 25% de l’aigua que s’extreu per regadiu. En el context del canvi climàtic amb canvis en els patrons de precipitacions i augment de la temperatura, les aigües subterrànies poden tenir un efecte directe sobre la disponibilitat d’aquest tipus de recurs a curt i llarg termini. De manera que la necessitat d’optimitzar-ne el seu ús i protegir-lo de la contaminació i la sobreexplotació és creixent. D’altra banda, l’ús cada cop més extensiu de noves tecnologies i la implementació de tècniques d’aprenentatge artificial així com l’augment de la capacitat computacional han facilitat la proliferació de projectes i tot tipus d’iniciatives que incorporen aquestes noves tècniques d’anàlisi i representació de dades al món de les aigües subterrànies (també anomenat hidrogeologia) tenint en compte la seva variabilitat espacial. En aquest sentit, l’ICGC (Institut Cartogràfic i Geològic de Catalunya) desenvolupa des de l’any 2010 el projecte de cartografia hidrogeològica de Catalunya a escala 1:25.000 que inclou la recollida sistemàtica de dades de camp sobre l’estat quantitatiu i qualitatiu de les aigües subterrànies en punts d’aigua (pous i fonts) i en concret anàlisis químiques de diverses tipologies per cadascun dels aqüífers cartografiats. Les dades de camp recollides i els resultats de les anàlisis químiques que es realitzen en el marc d’aquest projecte es recullen, validen i introdueixen a la Base de Dades del Mapa Hidrogeològic de Catalunya (BDSIMHCat). Institut Cartogràfic i Geològic de Catalunya (icgc.cat) que és l’eina de gestió de la informació hidrogeològica generada en els projectes d’hidrogeologia i geotèrmia de l’ICGC. La base està dissenyada i implementada en codi lliure mitjançant el sistema gestor de bases de dades relacional PostgreSQL juntament amb el mòdul PostGIS per afegir capacitats espacials. Les dades estan estructurades de manera que permet emmagatzemar-les, preservar-les i consultar-les en el marc dels projectes hidrogeològics desenvolupats. Així doncs es fa necessari una eina visual i intuïtiva per tal que usuaris externs puguin explotar aquestes dades per a poder treure resultats per als seus informes o plantejar noves hipòtesis, etc. Aquest treball s’estructura de la manera següent: primerament, es detallarà la metodologia que s’ha seguit per a desenvolupar l’aplicació, es donarà una descripció dels paquets R utilitzats, seguidament s’explicarà cada apartat i element que apareix a l’aplicació, i finalment les eines i transformacions estadístiques emprades. Després es presentaran els resultats d’un cas pràctic que és d’especial interès perquè fent servir l’aplicació s’obtindrà uns resultats molt similars a un estudi previ realitzat per l’Institut Cartogràfic. Per últim, es presentaran les conclusions finals i una discussió. 1 2 Objectius L’objectiu principal d’aquest treball és dissenyar una eina o interfície gràfica interactiva per tal de visualitzar la Base de dades del Mapa Hidrogeològic de Catalunya(BdMHCAt), mitjançant mapes i gràfics adequats. Com a objectiu secundari, fer ús de tècniques escaients per a analitzar el tipus de dades específiques com són les característiques dels pous i fonts d’aquest projecte. Finalment, un segon objectiu secundari serà realitzar els passos previs necessaris de processament de les dades (control de qualitat, transformacions, etc.). 2 3.2.4 Anàlisi univariant A la pàgina “univariant” apareixen els diferents violinplots (una combinació de boxplot i gràfic de densitat) de les variables pH, Redox, conductivitat elèctrica i la temperatura de l’aigua dels punts seleccionats. Al costat tenim els violinplots dels ions. A la figura 6 es mostra el contingut de l’“Anàlisi univariant” on apareix la pestanya “classes”. En aquesta pàgina es representen els ions estratificats per la variable “classes” (grups de pous o fonts trobats a partir d’un anàlisi de mescles gaussianes que s’explicarà més endavant). Els ions també queden estratificats utilitzant les seves dades originals o la seva transformació ILR. Aquesta transformació s’explica en més detall en els apartats següents. Figura 6: Pàgina univariant. 3.2.5 Anàlisi bivariant A la pàgina bivariant es mostren les relacions bivariants dels punts seleccionats. Els punts estan estratificats segons l’acrònim del seu aqüífer. En un principi la pàgina està buida, però a la barra lateral es poden seleccionar els gràfics que l’usuari vol veure. És important que, pel correcte funcionament de l’aplicació, aquesta pàgina inicialment quedi buida i que sigui l’usuari el que decideixi quins gràfics vol que apareguin. Si no es fes així, l’aplicació tardaria molt a carregar o, fins i tot, deixaria de funcionar a causa de la quantitat de gràfics que hauria de calcular. A la figura 7 es veu un exemple on se seleccionen els gràfics bivariants que es veuran en l’aplicació. En aquest cas la quantitat de calci contra la quantitat de sodi i la quantitat de calci contra la quantitat de magnesi dels punts seleccionats i estratificats per l’acrònim de l’aqüífer. 9 Figura 7: Pàgina bivariant. 3.2.6 Anàlisi multivariant En la pàgina “Anàlisi multivariant” de l’aplicatiu es poden seleccionar de forma interactiva els ions que s’utilitzaran per construir el model multivariant i, també, el nombre de clústers màxim que hi pot haver, tal com es veu en la figura 8. Si es prem “actualitza”, s’obté el millor model de mescles gaussianes, d’acord amb el criteri BIC, dels punts seleccionats. Sobre aquest model s’aplica una reducció de la dimensió en les dues primeres components principals (Dir1, Dir2). Així aconseguim el diagrama de dispersió dels punts seleccionats amb la seva respectiva classe o grup classificat. Per acabar, es pot descarregar la base de dades on hi ha la classe dels punts seleccionats com una nova columna al final. Figura 8: Pàgina multivariant. 10 A la figura 8 es mostren els dos diagrames de dispersió: el de l’esquerra s’ha tingut en compte que les dades són de tipus composicional i s’ha fet una transformació CLR (que és una transformació pròpia de les dades composicionals i que s’explicarà a continuació) a les dades prèviament, mentre que en el de la dreta s’ha fet amb les dades estandarditzades. D’aquesta manera es pot veure què succeeix si s’analitzen les dades sense tenir en compte que són de tipus composicional. En aquest cas, el diagrama de la dreta dona classificacions o grups espuris. A banda de resultats addicionals que proporciona l’anàlisi de mescles gaussianes, es poden veure les distribucions normals segons les dos primeres components (Dir1, Dir2) estratificades per les classes obtingudes, més informació del model i els gràfics BIC de tots els models. 3.2.7 Connexió entre les pàgines de l’aplicació Una de les funcionalitats més importants de l’aplicació és la següent: Si es clica sobre la frase “Actualitza BBDD amb el resultat del clustering” tots els resultats de totes les altres pàgines s’actualitzen d’acord amb les classes o grups trobats a partir del model de mescles gaussianes de la pàgina d’“Anàlisi multivariant”, tal com es veu a la figura 9. Ara els punts seleccionats del mapa estan estratificats segons la classe. Aquesta actualització afecta a tota l’aplicació: tots els gràfics estaran estratificats d’acord amb la variable classe i a la base de dades s’incorpora també aquesta variable. Figura 9: Pàgina principal mapa actualitzat. 11 3.3 Eines estadístiques 3.3.1 Dades composicionals Les dades composicionals (Peter Filzmoser 2009) són aquelles dades que representen una composició d’un tot. En el nostre cas, les variables hidroquímiques Ca, Mg, Na, K, Cl, CO3, SO4 i NO3, que utilitzem en el model de mescles gaussianes, són de tipus composicional perquè ens indiquen la composició química de la mostra. El punt diferencial, a part que les variables no són independents, és que la seva suma és una quantitat fixa no podem aplicar els mètodes estadístics habituals. S’ha d’aplicar diferents transformacions sobre aquestes variables. Si no es fes així, es podria obtenir relacions espúries i podria portar a conclusions incorrectes. La definició matemàtica formal de l’espai mostral d’unes dades composicionals és la següent: SD={X= [x1, ..., xD]∈RD|xi>0, i = 1, ..., D; D X i=1 xi=k} On D és el nombre de components, X és el vector positiu dels components xi, és la dada composicional i, i finalment k és la constant del sumatori de tots els components. Centered log ratio transform La transformació CLR resol el problema de dependència de les variables. Utilitzant la següent fórmula s’aconsegueix transformar les dades en independents i no correlacionades. En aquesta transformació es fa servir un logaritme i una operació de centrat que en aquest cas és la mitjana geomètrica. clr(x) = ln x1 n qQD i=1 xi , .., ln xD n qQD i=1 xi Finalment, aquest nou espai de dades logarítmic és simètric i la seva distribució és normal, la qual cosa permet fer ús de l’algorisme de mescles gaussianes. Isometric log-ratio transform Per a l’anàlisi univariant es fa servir aquesta transformació. Tant la tranformació ILR com en la transformació CLR persegueixen el mateix objectiu: tractar les dades composicionals. Però en el cas univariant és preferible l’ILR perquè conserva l’estructura composicional de les dades, la qual cosa permet fer comparacions relatives. En canvi, el CLR en estar en una formulació logarítmica perd aquestes propietats. La tranformació ILR té la següent fórmula: ilr(x) = r1 2ln(x 1−x) 3.3.2 Model de mescles gaussianes Un cop fet la transformació CLR sobre les dades, ja es pot aplicar l’anàlisi de mescles gaussianes (Fraley and Raftery 2002). Es tracta d’un algorisme no supervisat que serveix per classificar les dades en diferents 12 clústers. El model assumeix que, donat o condicionat a un grup, les dades o característiques d’un individu segueix una distribució normal multivariant, i l’objectiu és classificar cada individu a un grup donades les seves característiques. Així, la funció de versemblança de la mostra és la següent: n Y i=1 G X k=1 τkϕk(xi|µk,Σk)(1) On x representa les dades, G és el nombre de components o classes, τkés la probabilitat que l’observació pertany al k component. Per últim, ϕk(x|µk,Σk)és la funció de densitat normal multivariant del component k. ϕk(x|µk,Σk) = (2π)− p 2|Σk|− 1 2exp(−1 2(xi−µk)TΣ−1 k(xi−µk)(2) On µkiΣksón respectivament la mitjana i la desviació típica de la classe o grup k. L’algorisme per a classificar els individus de la mostra als diferents grups següeix els següents passos: Pas 1 Fixar el nombre de clústers (k) Pas 2 Inicialització dels centroides, k punts, de forma aleatòria Pas 3 Estimació dels paràmetres (mitjanes µi variàncies σ) amb el criteri de màxima versemblança (MV). Cal notar que la funció de versemblança de la mescla de gaussianes conté una suma dins d’un productori i, per tant, no té una solució tancada i calen mètodes iteratius com ara l’EM que té implementat el paquet mclust. Pas 4 Assignar cada individu de la mostra al clúster més probable seguint la regla de Bayes. D’aquesta manera s’assigna al grup amb la probabilitat a posteriori més gran, on aquesta probabilitat es calcula com: P(grupk|x) = P(x|grupk)P(grupk) PG j=1 P(x|grupj)P(grupj)(3) On x són els valors de les variables de l’individu (en el nostre cas serien les variables Ca, Mg, Na, K, Cl, CO3, SO4 i NO3) transformades amb CLR. Pas 5 Recalcular els centroides i els paràmetres utilitzant els punts assignats en l’apartat anterior Repetir 4i5fins a obtenir una convergència. Finalment, s’aplica una reducció de la dimensionalitat en dues components per poder visualitzar el model de mescles gaussianes estimat amb la seva respectiva classe o grup. Tot aquest procés es pot aplicar si se sap del nombre de classes o grups. Si no se sap, es pot efectuar amb diferent nombre de grups (per exemple de 2 a 6 grups) i escollir aquell que proporciona un BIC millor. 13 4 Cas pràctic Tot seguit es presentarà un cas particular que és d’especial interés ja que s’intentarà replicar el que es va fer en un estudi previ realitzat per l’Institut Cartogràfic (Camps Clemente Victor 2021). Primer de tot se seleccionen els polígons de diferents aqüífers de basament que inclouen tant punts d’aigua interiors (terra endins) com punts d’aigua més propers a la línia de costa i, per tant, susceptibles de ser afectats per intrusió salina (flux subterrani d’aigua salada provinent del mar), en un àmbit situat a l‘Alt Empordà (NE de Catalunya). Figura 10: Pàgina principal mapa cas pràctic. En l’anàlisi univariant (gràfic de la dreta), s’observa una gran dispersió en el contingut de clorur (Cl) i Sodi (Na) indicatius del nivell de salinitat de l’aigua. Aquesta dispersió probablement indica que, de les mostres seleccionades, algunes probablement estan afectades per intrusió salina o aerosol d’aigua marina. Figura 11: Pàgina univariant cas pràctic. 14 Quan passem a l’apartat d’anàlisi multivariant amb el model de mescles gaussianes, el biplot obtingut de les dues primeres components principals aplicant les transformacions CLR sobre els continguts dels ions majoritaris en l’aigua (calci, sodi, potassi, bicarbonats, clorur, sulfats i nitrats), es veuen tres grups clarament diferenciats (1, 2 i 3) (Figura 12, imatge de l’esquerra). En canvi, el biplot de coordenades on les dades estan normalitzades i estandarditzades (log z-score) (sense aplicar la transformació CLR), es veuen dos grups, però que no estan gaire diferenciats i força solapats. Figura 12: Pàgina multivariant cas pràctic. A la figura 13 es representen els punts en la primera coordenada en forma de gràfics de densitat estratificats segons la classe. Aquesta diferència entre els resultats de les dades aplicant les transformacions CLR i els resultats aplicant les transformacions log z-score, és encara més evident. S’observa bé la separació entre els tres clústers pel model CLR (a l’esquerre) respecte als dos clústers del model log z-score (aplicar un logaritme a les dades originals i després restar la mitjana i dividir el resultat per la desviació típica) on se superposen molt (a la dreta). Figura 13: Pàgina multivariant cas pràctic. 15 A la figura 14 es mostren els punts seleccionats sobre el mapa i en el diagrama de Piper diferenciats segons els clústers obtinguts a partir del model de mescles gaussianes amb la transformació CLR. S’observa com les mostres classificades dins el clúster 1 (de color vermell) són aquelles amb major contingut de clorurs i sodi (ions Cli Na+) i que aquestes se situen properes al mar. Concretament, i segons el diagrama de Piper, les mostres del clúster 1 es correspondrien a mostres cloururades-sòdiques i, per tant, molt possiblement afectades per intrusió salina. Figura 14: Pàgina principal mapa actualitzat cas pràctic. Cal tenir en compte que la intrusió salina pot ser un procés natural o bé afavorit antròpicament pel bombament de pous d’aigua dolça propers a la línia de costa i que, per tant, per determinar l’origen de la salinització dels pous, caldria realitzar estudis de més detall sobre la dinàmica dels fluxos d’aigua subterranis i l’existència d’extraccions properes. En qualsevol cas, el resultat final obtingut és molt similar i coherent amb els resultats d’un estudi fet per l’Institut Cartogràfic i Geològic de Catalunya en el marc del projecte GeoERA TACTIC anomenat “Assessment of salt-/sea water intrusion status and vulnerability”. En aquest projecte (Camps Clemente Victor 2021) es va determinar el llindar pel que fa a la concentració de clorurs a partir del qual una zona es pot considerar afectada per intrusió salina i s’avaluà l’evolució històrica d’aquesta àrea afectada en l’àmbit de l’Alt Empordà. Com es pot veure a la següent imatge, on es delimita l’aqüífer fluviodeltaic del Fluvià i la Muga, els àmbits amb major concentració de clorurs en l’aigua subterrània en els anys 1995, 2011 i 2017 coincideix significativament amb la localització dels punts que l’aplicació desenvolupada en el seu apartat d’anàlisi multivariant s’han classificat segons el clúster 1 tenint en compte la transformació de les dades amb el mètode CLR. Terra endins, les mostres classificades dins els clústers 2 i 3, es corresponen majoritàriament a aigües bicarbonatades càlciques amb un contingut baix de clorur i sodi coincidint amb la zona no afectada per intrusió salina. 16 17 5 Conclusions L’aplicació final ha estat publicada al servidor que ofereix shiny (shinyapps.io), i que permet als usuaris de R publicar i compartir les aplicacions web desenvolupades utilitzant R shiny. En el següent enllaç hi ha l’aplicació final. https://icgcshiny01022022.shinyapps.io/shinyICGC/. Els objectius principals d’aquest treball queden assolits ja que s’ha aconseguit crear una aplicació interactiva per facilitar la consulta i l’anàlisi estadístic de les dades de l’Institut Cartogràfic i Geològic de Catalunya, tenint especialment en compte que algunes d’aquestes variables són de tipus composicional. Aquesta aplicació va destinada a les necessitats de l’Institut Cartogràfic on l’usuari no necessita experiència en R, però sí tenir coneixements geològics per poder interpretar bé els resultats i s’ha fet a mida i demanda dels responsables de l’Institut Cartogràfic. Com a conclusió o reflexió final, no només és necessari que els aplicatius Shiny facin servir les eines estadístiques adients, cosa que és relativament fàcil perquè R ofereix una quantitat molt gran de funcionalitats en aquest sentit, sinó que l’aspecte visual de l’aplicatiu web ha de ser atractiu per l’usuari i ha de concordar amb la imatge que vol projectar la institució o l’empresa. Això només es pot aconseguir fent servir codi CSS i HTML en l’aplicació. El procés de realització d’aquell treball m’ha servit per aprendre o acabar d’aprendre com es confecciona una aplicació web amb Shiny, crear i personalitzar mapes i també conceptes propis d’aquest àmbit o camp com són les coordenades i les seves transformacions o les capes, etc., diagrames de Piper, tractament de dades composicionals i aprofundir amb models de classificació no supervisats com són les mescles gaussianes. La major part del volum de feina d’aquest treball ha estat la confecció de l’aplicació amb les eines de Shiny i que no queda del tot reflectit en aquesta memòria, on presento totes les funcionalitats i les eines estadístiques utilitzades. Les eines emprades per fer aquesta aplicació no s’ensenyen en el màster com, per exemple, la creació del mapa fent servir leaflet, o la confecció d’aplicatius amb Shiny. Com a treball futur quedaria pendent perfeccionar i afegir noves funcionalitats a l’aplicació actual com ara la possibilitat que l’usuari pugui analitzar la seva pròpia base de dades de pous i fonts. També seria interessant provar amb més exemples de casos pràctics, i confeccionar un extens manual d’ajuda. 18