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

1. Introducción

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.

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:

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

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.

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 .

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.

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.

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.

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 km².

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

8. Estilizar la capa del MDE

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

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

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.

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

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.

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 .

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.

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.

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.