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, 9 de mayo de 2025, 09:14 |
Descripción
Tabla de contenidos
- 1. Introducción
- 2. Cargar datos de pozos desde un GeoNode IDE (Infraestructura de Datos Espaciales)
- 3. Exportar datos de la capa WFS hacia el archivo GeoPackage
- 4. Seleccionar pozos de un área específica
- 5. Descargar el MDE del SRTM 1-Arc Second
- 6. Calcular la densidad de pozos
- 7. Cortar y reproyectar el Modelo Dígital de Elevación (MDE)
- 8. Estilizar la capa del MDE
- 9. Muestrear las elevaciones al MDE a la ubicación de los pozos
- 10. Comparar los datos de elevación reportada con los valores del MDE
- 11. Reemplazar valores NoData con los valores de elevación MDE
- 12. Calcular el nivel del agua subterránea en los pozos
- 13. Interpolar los niveles de agua subterránea de los pozos a un ráster
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.
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
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.
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 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...
8. Estilizar la capa del MDE
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 .
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.
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.