Tutorial: Creación de un mapa del nivel de agua subterránea a partir de los datos de pozos y un Modelo Dígital de Elevación (DEM).

Sitio: OpenCourseWare for GIS
Curso: Capacitación SIG (Sistema de Información Geográfica) orientado a aplicaciones hidrogeológicas.
Libro: Tutorial: Creación de un mapa del nivel de agua subterránea a partir de los datos de pozos y un Modelo Dígital de Elevación (DEM).
Imprimido por: Guest user
Día: viernes, 29 de marzo de 2024, 09:50

Descripción

1. Introducción

Al terminar este tutorial será capaz de:

  • Cargar datos de pozos desde una IDE (Infraestructura de Datos Espaciales) en QGIS
  • Calcular la densidad de los pozos dentro de un polígono
  • Descargar un Modelo Dígital de Elevación (DEM)  SRTM 1-Arc Second
  • Cortar y reproyectarClip and reproject a DEM
  • Estilizar un DEM
  • Muestrear la elevación de un pozo y compararlos con el atributo de elevación
  • Correguir cálculos en la tabla de atributos
  • Interpolar los niveles de agua subterránea de los pozos a un ráster

2. Cargar datos de pozos desde un GeoNode IDE (Infraestructura de Datos Espaciales)

Iniciaremos con la carga de datos de pozos desde un GeoNodo IDE a QGIS.

En este tutorial utilizaremos la base de datos Borehole Database in the Stampriet Transboundary Aquifer desde Orange-Senqu River Basin GIS Server, el cual es el GeoNode SDI.

1. Verifica el metadato de Borehole Database in the Stampriet Transboundary Aquifer.  Selecciona la pestaña de Atributtes y selecciona los atributos en la capa.

stasboreholes in GeoNode

2. Ejecuta  QGIS con un proyecto vacío.

3. Da un click en el botón Open Data Source Manager   en la barra de herramientas.

4. Selecciona la pestaña GeoNode.

5. Da un click en botón New para crear un nuevo servicio de conexión.

6. En la ventana de diálogo Create a New GeoNode Connection teclea en el campo de nombre ORASECOM y copia el URL http://gis.orasecom.org.

7. Da un click en el botón Test Connection.

Si la prueba es éxitosa observarás está ventana:


Si la conexión falla, verifica la conexión a Internet y el URL.

8. Da click en el botón OK para cerrar la ventana popup.

9. Da click en el botón OK en la ventana de diálogo  y una nueva conexión será agregada.

10. Da un click en el botón Connect .

11. Sí aparece un error en una ventana puedes ignorarla. Da click en el botón OK para borrar el mensaje de error.

Observarás las capas en lista del GeoNodo.

Utilizaremos los datos del Acuífero Stampriet.

12. Teclea stampriet como filter.

13. Da un click en la capa stasboreholes WFS (mostrada en la figura inferior).

Utilizaremos el servicio WFS (Web Feature Service) , porque requerimos acceder a los datos y no sólo a una imagen de los mismos.

14. Da un click en el botón Add  que agrega la capa al  "map canvas".

15. Da un click al botón Close para cerrar la ventana de diálogo.


Tardará un poco en cargarse desde Internet.


Espere hasta que cargue la capa completamente. Entonces observarás esto:

Loading the stasboreholes


En la siguiente sección vamos a guardar esta capa WFS en un GeoPackage para utilizarlo localmente.


3. Exportar datos de la capa WFS hacia el archivo GeoPackage

Nuestro conjunto de datos aún está en formato WFS en un Sistema Coordenado Geográfico con coordenadas de latitud/longitud en grados. Para nuestro propósito de interpolación de la elevación, es necesario, almacenar el archivo localmente y reproyectarlo en UTM Zone 34S / WGS-84.


1. En el panel de Layers da un click derecho en la capa geonode:stasboreholes .

2. Selecciona Export | Save features as...

3. En la ventana de diálogo Save Vector Layer as... selecciona GeoPackage como tipo de formato Format. Da un click en el botón  para navegar en  los archivos del disco duro donde se guardará el archivo GeoPackage (e.g. Z:\Stampriet). Nombra al archivo como Stampriet_Data.gpkg.  Teclea el nombre de la capa  (Layer name ) como Boreholes.

4. Da in click en el botón  que cambie la proyección de la capa de salida.

5. En la ventana de diálogo  Coordinate Reference System Selector selecciona 32734 (que es el código EPSG ) en el campo del filtro (Filter)  en la búsqueda de la proyección WGS 84 / UTM Zone 34S . Da un click sobre el nombre de la proyección y da un click en el botón  OK  para regresar a la ventana de diálogo.

6. Da un click en el botón OK que ejecuta la exportación.

7. Después la capa Stampriet_Data Boreholes  se agregará al "Map Canvas" , borra la capa geonode:stasboreholes  haciendo un click derecho sobre esta y selecciona Remove layer.... Confirma con el botón OK.


Ahora tendremos una copia local de la capa de pozos.

8. Cambia la proyección OTF del proyecto: en la panel  de "Layers" panel  da un click derecho sobre la capa Stampriet_Data Boreholes  y selecciona Set CRS | Set Project CRS from Layer

9. Salva el proyecto.

En la siguiente sección, seleccionaremos los pozos de un área en específico.


4. Seleccionar pozos de un área específica

En este tutorial, sólo nos interesan los pozos en el área de estudio alrededos de Stampriet.

Un polígono imaginario es proporcionado en los datos del tutorial.

1. Agrega el shapefile  Stampriet_boundary  al proyecto.

2. En el panel Layers da un click derecho sobre la capa Stampriet_boundary y elige la opción Zoom to layer.


3. Estiliza la capa con un contorno rojo.

Stampriet boundary styled

Ahora, cortaremos  de la capa de pozos, aquellos que se encuentren dentro de los límites y guardar los resultados al archivo GeoPackage.

4. Desde el menú principal elige Vector | Clip...

5.  En la ventana de  diálogoIClip selecciona Stampriet_Data Boreholes como una capa de entrada (Input layer ) y la capa Stampriet_boundary como capa sobrepuesta (Overlay layer). Guarda el resultado del geoproceso Clipped al archivo Geopackage Stampriet_Data con el nombre de la capa boreholes_clipped.

6. Da un click en el botón Run. Cierra la ventana de diálogo  después del geoprocesamiento, da un click en el botón Close.

7. Borra la capa Stampriet_Data Boreholes desde el panel Layers .

Selected boreholes Stampriet


En la siguiente sección, calcularemos la densidad de pozos en está área.


5. Descargar el MDE del SRTM 1-Arc Second

En esta sección, agregaremos un MDE al proyecto, Utilizaremos el producto Shuttle Radar Topography Mission (SRTM) 1 Arc Second. Este es un MDE de acceso libre que cubres 56°S a 60°N y tiene una resolución espacial de aproximadamente 30 m sobre el Ecuador.

Puedes descargar los tiles de SRTM  desde el sitio USGS Earth Explorer y/o a través del complemento. Emplearemos el complemento "SRTM Downloader" en QGIS.

1. Desde el  menú principal, selecciona Plugins | Manage and Install Plugins...

2. Busca e instala el complemento SRTM-Downloader .

3. En el menú principal, da un click en el icono SRTM-Downloader   que ejecuta el complemento.

4. En la ventana de diálogo SRTM-Downloader da un click al botón Set canvas extent (asegurate que tienes la extensión del área de estudio!).  La ruta del archivo de salida, navega al folder donde se encuentran los datos del proyecto.

5. Da un click en el botón Download.

6. La primera ejecución, abrirá una ventana que solicitá la información de acceso Login.  Se requiere una cuenta en https://urs.earthdata.nasa.gov/users/new que permite la descarga de los Tiles SRTM.

7. Introduce Username y Password. Verifica el checkbox Save Credentials  y da click en el botón OK.

Ahora, inicia la descarga de los Tiles SRTM que cubren el área de estudio.

8. Da un click en el botón OK  en el popup.

9. Cierra la ventana de diálogo.

Como puedes observar, el Tile SRTM es mayor que el área de estudio. En la siguiente sección , cortaremos y reproyectaremos el MDE.

SRTM tile added to project


6. Calcular la densidad de pozos

Ya tenemos los pozos dentro del área de estudio (o bien dentro de un acuífero) donde calcularemos la densidad de pozos.

1. Abre la caja de procesamiento (Processing Toolbox)  desde el menú principal, elige Processing | Toolbox.

2. En la caja de procesamiento (Processing Toolbox) selecciona Vector analysis | Count points in polygon.

3. En la ventana de diálogo Count Points in Polygon selecciona la capa Stampriet_boundary en la sección de Polygons y la capa boreholes_clipped en la sección de Points. Cambia en el campo Count field name a No of boreholes. Guarda la salida al archivo   GeoPackage Count  Stampriet_Data  con el nombre de la capa Stampriet study area.

Usamos ese nombre de archivo de salida, porque la herramienta hace un copia del polígono Stampriet_boundary y agrega el campo "No of boreholes". Aún no tenemos la capa del área de estudio en el archivo GeoPackage, así que lo agregaremos  de la siguiente forma.

4. Da click en el botón Run. Cierra (Close ) la ventana después del geoprocesamiento.

Considera que si la capa de entrada tiene múltiples polígonos (e.g. diferentes acuíferos) producirá el conteo por polígono. Puedes seguir los pasos para calcular la densidad de pozos en todos los polígonos.

5. Copia el estilo de la capa Stampriet_boundary a la capa Stampriet study area: en el panel Layers da click derecho sobre la capa Stampriet_boundary y selecciona Styles | Copy Style | All Style Categories.  Entonces da un click derecho sobre la capa Stampriet study area  y elige Styles | Paste Style | All Style Categories.

6. Borra la capa Stampriet_boundary del panel Layers.

Ahora, revisemos la tabla de estudio de la capa Stampriet study area .

7. En el panel Layers da un click derecho sobre la capa Stampriet study area y selecciona Open Attribute Table.

Observarás que hay 64 pozos en el área de estudio.

Ahora, calcularemos el área de la zona de estudio.

8. En la tabla de atributos da un click en el botón  que active la edición.

9. Da un click en el icono  que abre la calculadora de campos (Field calculator).

10. En la ventana de la calculadora de campos (Field Calculator)  crea un nuevo campo (Create a new field) con el nombre del campo de salida (Output field name) como Area (km2) y  el tipo de dato del campo de salida (Output field type) Decimal number (real).

11. Bajo la sección de la Expression (expresión) escoge la función en el panel medio  de la parte Geometry | $area da un doble click sobre $area que lo agrega a la sección de la expresión, La variable $area  calcula el área del polígono en unidades del mapa.  En nuestro caso en metros cuadrados. Completa la expresión  de la siguiente forma:
$area / 1000000

Utiliza el botón  que indica la división.

12. Da un click en el botón OK.

El resultado lucirá como la figura inferior.

El siguiente paso, es el cálculo de la densidad de pozos.

13. Da un click en el icono  que abre la calculadora de campos (Field calculator).

14. En la ventana de diálogo de la calculadora de campos (Field Calculator ), verifica el checkbox Create a new field , en el campo Output field name  utiliza el nombre Borehole density per km2  y el campo de salida (Output field type) de tipo Decimal number (real). En la sección de la expresión ingresa :

 "No of boreholes"  /  "Area (km2)"

Recuerda que puedes agregar los campos, haciendo doble click en el nombre de los campos bajo la sección  Fields and values en el medio del panel. Usa el botón  que agrega el símbolo de división.

15. Da un click en el botón OK.

16.  Habilitá la edición, activando el icono   y guarda con el botón Save.

Calculemos la densidad de pozos en el área de estudio que es 0.2 per km2.


En la siguiente sección, agregaremos un Modelo Dígital de Elevación (MDE) al proyecto.


7. Cortar y reproyectar el Modelo Dígital de Elevación (MDE)

El tile del SRTM es agregado al  Map Canvas. Si te acercas a la capa S25E018.hgt layer observarás que cubre un área amplia.  Asimismo el Sistema Coordenado Geográfico (GCS) con coordenadas de  latitude/longitude en grados (EPSG:4326).

Así, recortaremos y reproyectaremos el MDE y se agregará al archivo GeoPackage.

1. En el panel  Layers haz un click derecho sobre la capa S25E018.hgt y selecciona la opción Export | Save As...

2.  En la ventana de diálogo  Save Raster Layer as... selecciona al archivo GeoPackage como formato de salida. Navega al archivo Stampriet_Data.gpkg y selecciona el DEM como nombre de la capa (Layer name).  Cambia el  CRS a Project CRS: EPSG:32734.
Bajo el sección  de Extent elige el botón Calculate from Layer y selecciona la capa Stampriet study area.
En la parte de Resolution cambia los valores de la resolución Horizontal y Vertical a 30 m.
Verifica el checkbox a  No data values y agrega el valor de -9999 bajo FromTo.  Este será el valor  en las celdas fuera de la zona de estudio, después del corte y reproyección.


3. Deja los demás parámetros de configuración como "default" y da click en el botón OK.

4. Borra la capa S25E018.hgt desde el panel Layers y acercate a la zona de estudio con los resultados.

En la siguiente sección, estilizaremos el MDE.

8. Estilizar la capa del MDE

Para una buena visualización estiliza el MDE con una rampa de color y un hillshading.

1. Da un click derecho sobre el MDE y selecciona Duplicate Layer.

2. Renombra la capa MDE copy  como Hillshade.

3. Selecciona la capa MDE  en el panel Layers y da un click en el icono  que abre el panel de estilizado (Layer Styling).

4. Cambia el tipo de renderer a  Singleband pseudocolor.

5. Da un click en la rampa de color (Color ramp) sobre la flecha y selecciona Create New Color Ramp...


6. En la ventana de diálogo Color ramp type escoge Catalog: cpt-city y da un click en el botónOK.


7. En la ventana de diálogo de Cpt-city Color Ramp elige una paleta de colores de la categoría Topography, e.g. elevation y da click en el botón OK.


8.  En el panel de estilizado de la capa (Layer Styling) selecciona la opción de Classify en la rampa de color.

9. Bajo la sección Layer Rendering cambia del modo Blending mode a Multiply.  Esto es necesario para el siguiente paso, donde diseñaremos el sombreado que veremos a tráves del MDE.


10. Desde el panel Layers activa la capa Hillshade .

11. En el panel  Layer Styling asegurate que se estiliza la capa Hillshade . Cambia el tipo de renderer a Hillshade. Cambia en la sección de remuestreo , los parámetros de ResamplingZoomed in  a  Bilinear y   Zoomed out a Average. Se creará una visualización más suave. Mantén las demás configuraciones por default.

Ahora, el proyecto luce como se muestra en la imagen inferior.

DEM styled

En la siguiente sección, agregaremos los valores de las elevaciones del MDE a la tabla de atributos de los pozos.


9. Muestrear las elevaciones al MDE a la ubicación de los pozos

Muestrearemos las elevaciones del MDE a la ubicación de los pozos.

Para ello, se requiere la instalación del complemento Point sampling tool . Considera que requieres una conexión a Internet que permita la instalación del complemento.

1. Desde el menu principal dirígete In the main menu go to Plugins | Manage and Install Plugins...

2. En la ventana de diálogo de  Plugins busca  Point sampling tool  e indica Install Plugin.

Install Point sampling tool

3. Da un click al botón  desde la barra de herramientas y ejecuta  el complemento Point Sampling Tool 

4.  En la ventana de diálogo Point Sampling Tool en la pesptaña General selecciona la capa boreholes_clipped como Layer containing sampling points. Selecciona todos los items con el nombre de boreholes_clipped y la capa DEM . Considera que estas capas deben habilitarse en el panel Layers  en la forma en que aparecerá aquí. Guarda la capa como  Output point vector layer a un archivo tipo ESRI Shapefile, porque esta herramienta no puede agregarse a un archivo GeoPackage. Nombra al archivo como Boreholes_Z.shp.

5. Ahora, selecciona la pestaña Fields .  Puedes cambiar el nombre de los campos de salida si es necesario,  sólo cambiaremos el campo del  Z del MDE.

6. Da click en el botón OK. Cierra la ventana despúes del geoprocesamiento.

7. Verifica la tabla de atributos de la capa  Boreholes_Z .

Verifica que el campo Z ha sido agregado.

En la siguiente sección, crearemos una gráfica scatter de la elevación y el valor Z desde el MDE.


10. Comparar los datos de elevación reportada con los valores del MDE

En la tabla de atributos, el campo Boreholes_Z  contiene el valor de la elevación del DEM. Comparemos los dos campos por medio de una gráfica scatter.

El campo de elevación contiene valores NoData. No consideraremos estos, en las gráficas scatter.

1. Abre la tabla de atributos de la capa Boreholes_Z.

2. Da un click en el botón Select by attributes  .

3. Escribe la expresión:

"elevation" < 0

Select elevation < 0

4. Da un click en el botón Select Features y después cierra con el botón Close la ventana.

10 entidades han sido seleccionados. Borraremos desde la tabla de atributos con la habilitación del modo de edición y despuésda un click en el icono "trash". En el caso de la gráfica scatter, invertaremos la selección y elegiremos sólo a las entidades seleccionadas.

5. Da un click en el botón   que invierte la selección.

Elaborar la gráfica scatter requerimos el complemento Data Plotly .

6. Instala el complemento Data Plotly desde el administrador de complementos (Plugins manager).

7. Da un click en el icono  desde la barra de herramienta que abre el panel Data Plotly .

8. En el panel Data Plotly selecciona Scatter Plot como  Plot type. Escoge Boreholes_Z.shp como  Layer. Verifica el checkbox  de  Use only selected features. Selecciona  elevation para el campo X field y  Z para el campo Y field. Mantén el estilizado por default.

9. Da un click en el icono   que es la sección de configuración de la gráfica.

10. Deshabilita el checkbox Show legend . Cambia el título de la gráfica Plot title a  Borehole vs DEM elevation. Borra el título de la leyenda Legend title. Escribe Borehole elevation (m) para la etiqueta de X X label y DEM elevation (m) para la etiqueta de Y label. Mantén el resto por default.

11. Da un click en Create Plot  en la parte inferior del panel.

Observa la gráfica.

Scatter plot elevations


En la siguiente sección reemplazaremos los valores faltantes de la elevación con los valores de la elevación del MDE.


11. Reemplazar valores NoData con los valores de elevación MDE

Debido a que el diagrama de dispersión ( Scatter plot) parece razonable, reemplazaremos los valores de nodata en el campo de elevación de la capa Boreholes_Z.shp con los valores del MDE.

1. Abre la tabla de atributos de la capa Boreholes_Z.shp.

2. Da click en el icono  que deselecciona atributos.

3. Da click al icono  que habilita el modo de edición.

4. Cambia la salida a la elevación y da un click en el icono  que permite construir la expresión.

5. En la ventana de diálogo de la expresión contruiremos la siguiente expresión:

if( "elevation" < 0,  "Z" ,  "elevation" )

Esto significa: Si el campo de elevación tiene un valor menor a cero, entonces utiliza los valores del campo Z, en caso contrario, utiliza el campo de elevación  (elevation).

6. Da click al botón OK.

La expresión debería lucir como esto en la tabla de atributos;


7. Da click en el botónUpdate All.

Ahora, todos los valores faltantes han sido reemplazados con los valores de elevación desde el MDE.

8. Da click en el icono   que deshabilita la edición y da click en el botón Save para guardar los resultados.


En la siguiente sección, calcularemos el nivel absoluto del agua subterránea en los pozos.


12. Calcular el nivel del agua subterránea en los pozos

En la tabla de atributos de los pozos tenemos la profundidad del manto freático.

En esta sección, calcularemos la elevación absoluta del manto freático de los pozos.

1. Abre la tabal de atributos de la capa Boreholes_Z.shp.

2. Da click al icono  que habilita el modo de edición.

Primero, borraremos las entidades con datos faltantes en el campo depth_m.

3. Da un click  en el icono que abre la ventana de  Select features using an expression .

4. Construye la siguiente expresión

"depth_m" < 0

5. Da un click en el botón Select Features. Cierra la ventana de diálogo.

6. Da un click en el icono  que borra los registros seleccionados.

7. Da un click en el icono   que guarda las ediciones.

8. Da un click en el icono   que abre la ventana de la calculadora de campos (Field calculator).

9. En la ventana de diálogo de la calculadora de campos (Field Calculator) crea un nuevo campo (Create a new field ) con el nombre del campo de salida  GW_level_m y  el tipo del campo de salida (Output field type) como Decimal number (real).  Mantén las otras configuraciones como default con la siguiente expresión:

"elevation" - "depth_m"

10. Da click en el botón OK.

11. Verifica los resultados y da un click al icono de edición   . Guarda los resultados con el botón Save .


En la siguiente sección, interpolaremos los niveles del agua subterránea a un archivo raster.

13. Interpolar los niveles de agua subterránea de los pozos a un ráster

Tenemos los niveles de agua subterránea en los pozos, podemos interpolar los puntos a un ráster.

Utilizaremos dos métodos: Inverse Distance Weighting (IDW)  y polígonos de  Thiessen.

Empecemos con la interpolación IDW.

1. En la caja de herramientas (Processing Toolbox) selecciona Interpolation | IDW interpolation.

2. En la ventana de diálogo IDW Interpolation selecciona la capa Boreholes_Z.shp como capa vectorial (Vector layer).  Usa el atributo GW_level_m como un atributo de interpolación (Interpolation attribute). Da click al icono  que agrega la configuración a la lista. Mantén el  parámetro Distance coefficient P a 2  para una función exponential . Selecciona el Extent de la capa Stampriet study area. Cambia el tamaño del pixel (Pixel size) a 30 m.  Guarda el ráster interpolado (Intepolated raster) a la carpeta del proyecto como  GW_level_IDW.tif.

3. Da click en el botón Run. Cierra la ventana después de que termite el geoprocesamiento.

Estiliza el resultado como una rampa de color.

4. Selecciona la capa GW_level_IDWen el panel Layers y da click en el icono  que abre el panel de estilización Layer Styling.

5. En el panel Layer Styling cambia al  renderer a Singleband pseudocolor  y elige All Color Ramps | RdBu

El panel Layer Styling luce así:

Y la capa interpolada luce como la figura inferior. Puedes jugar con diferentes rampas y clases que mejoren la visualización. El interpolador IDW crea una salida suave, porque la cantidad de pozos es menor, hay mayor sesgo hacia valores individuales y valores atípicos.


IDW interpolation result


Ahora, apliquemos el método de polígonos de Thiessen (aka nearest neighbour oVoronoi tesselation).

6. Desde el menú principal, selecciona Raster | Analysis | Grid (Nearest Neighbor).

7. En la ventana de diálogo Grid (Nearest Neighbor) selecciona la capa Boreholes_Z.shp como capa de puntos (Point layer). Bajo la sección  Advanced selecciona como valor Z value from field al campo  GW_level_m . Guarda la  salida de la interpolación (Interpolated (Nearest Neighbor)  al archivo GW_level_Thiessen.tif. Mantén el resto de las configuraciones como default.

8. Da un click al botón Run. Cierra la ventana después del geoprocesamiento.

9. Estiliza el resultado.

Este algoritmo no cubre el área completa, es limitada al método convex hull de la capa de puntos de entrada. El resultado es menos suave que la interpolación IDW. Sin embargo, con datos escasos, el valor más cercano podría ser la mejor aproximación.

Thiessen interpolation result

Hemos terminado el proyecto con la estilización de los pozos y la exportación de todas las capas que deseemos conservar al archivo  GeoPackage.