Crear GeoTiff de una cuenca con herramientas libres


Howto generate drainage basins from ASTER GDEM Project data.

Existen varias fuentes para obtener DEM de forma gratuita, pero la página de ASTER GDEM es la que me ha dado mejores resultados ya que entrega una precisión de 30 m. para Chile. Existe información adicional que se puede descargar desde WIST, donde podemos usar la misma cuenta creada para ASTER GDEM.

Lamentablemente para crear el shape de la cuenca aún no he podido investigar si se puede hacer con los módulos de Grass, como por ejemplo r.watershed , r.water.outlet , etc. Por lo que tuve que recurrir al software privativo WMS 7.1 , siguiendo el manual que nos pasaron en el curso de Hidrología CI41C. Este software trabaja archivo en formato GridFloat, por lo que aprovechamos bajar los DEM desde la USGS los cuales tienen una precisión de sólo 90 m. Con esto podemos crear 3 shapes: _arc.shp (líneas con los cauces) , _poly.shp (polígono con la divisoria de aguas) y _pt.shp de la cuenca (punto salida de la cuenca).

Con el polígono con la divisoria de aguas es posible “cortar” desde el Geotiff de la zona la información referente sólo a la cuenca. Para cual debemos seguir los siguientes pasos.

  • Descargar los geotiff de la zona desde la página de ASTER GDEM , los cuales descargaremos en pedazos de 1 grado x 1 grado, por ejemplo ASTGTM_S29W071_dem.tif contiene la información del cuadrante comprendido entre S29 – S28 , W071 – W072. Junto a los archivos _dem se encuentra un archivo _num que corresponde al error de la medición.
  • Ahora nos toca unir todos los pedazos en uno que incluya toda la cuenca. Para eso no ayudaremos con el comando gdal_merge

# Con el siguiente comando unimos los archivos ASTGTM_S29W068_dem.tif , ASTGTM_S30W071_dem.tif y ASTGTM_S32W069_dem.tif , etc. y los guardamos en el archivo elqui-mosaico.tif ,  en este caso los archivos .tif están en el mismo directorio donde estamos trabajando.
$ gdal_merge.py -o elqui-mosaico.tif ASTGTM_S29W068_dem.tif ASTGTM_S30W071_dem.tif ASTGTM_S32W069_dem.tif
# Para ver la información de los archivos .tif , en este caso gdalinfo usamos el archivo que generamos.
$ gdalinfo elqui-mosaico.tif

Driver: GTiff/GeoTIFF
Files: elqui-mosaico.tif
Size is 14401, 14401
Coordinate System is:
GEOGCS[“WGS 84”,
DATUM[“WGS_1984”,
SPHEROID[“WGS 84”,6378137,298.2572235630016,
AUTHORITY[“EPSG”,”7030″]],
AUTHORITY[“EPSG”,”6326″]],
PRIMEM[“Greenwich”,0],
UNIT[“degree”,0.0174532925199433],
AUTHORITY[“EPSG”,”4326″]]
Origin = (-71.000138888888884,-27.999861111111112)
Pixel Size = (0.000277777777778,-0.000277777777778)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( -71.0001389, -27.9998611) ( 71d 0’0.50″W, 27d59’59.50″S)
Lower Left  ( -71.0001389, -32.0001389) ( 71d 0’0.50″W, 32d 0’0.50″S)
Upper Right ( -66.9998611, -27.9998611) ( 66d59’59.50″W, 27d59’59.50″S)
Lower Right ( -66.9998611, -32.0001389) ( 66d59’59.50″W, 32d 0’0.50″S)
Center      ( -69.0000000, -30.0000000) ( 69d 0’0.00″W, 30d 0’0.00″S)
Band 1 Block=14401×1 Type=Int16, ColorInterp=Gray

Donde podemos ver mucha información interesante, como por ejemplo: Formato GTiff/GeoTIFF, Coordenadas geodésicas Datum WGS 84 , el origen, las coordenadas de las esquinas y más abajo las Bandas de información en este caso sólo 1 del tipo Int16, la cual correspondería a las elevaciones.
Podemos abrir nuestro mosaico verlo en Qgis, en propiedades cambiamos el mapa de color a pseudocolor y podremos ver algo como esto:

Imagen GeoTiff del sector

Ya tenemos nuestro archivo mosaico listo para recortar la información de la cuenca, momento …. no es tan fácil, resulta que cuando trabajamos en WMS transformamos nuestro DEM a WGS 84 / UTM zone 19S y los shapes de las cuencas quedaron en esa proyección.

# Para ver la información de los shapes, usamos el comando ogrinfo, en este caso la divisoria de aguas de la cuenca en el archivo estacion01_poly.sh:

$ ogrinfo -al -geom=NO estacion01_poly.shp

INFO: Open of `estacion01_poly.shp’
using driver `ESRI Shapefile’ successful.
Layer name: estacion01_poly
Geometry: Polygon
Feature Count: 4
Extent: (392889.046537, 6679377.052909) – (415780.512166, 6711932.273227)
Layer SRS WKT:
PROJCS[“WGS 84 / UTM zone 19S”,
GEOGCS[“WGS 84”,
DATUM[“WGS_1984”,
SPHEROID[“WGS 84”,6378137,298.257223563,
AUTHORITY[“EPSG”,”7030″]],
AUTHORITY[“EPSG”,”6326″]],
PRIMEM[“Greenwich”,0,
AUTHORITY[“EPSG”,”8901″]],
UNIT[“degree”,0.01745329251994328,
AUTHORITY[“EPSG”,”9122″]],
AUTHORITY[“EPSG”,”4326″]],
UNIT[“metre”,1,
AUTHORITY[“EPSG”,”9001″]],
PROJECTION[“Transverse_Mercator”],
PARAMETER[“latitude_of_origin”,0],
PARAMETER[“central_meridian”,-69],
PARAMETER[“scale_factor”,0.9996],
PARAMETER[“false_easting”,500000],
PARAMETER[“false_northing”,10000000],
AUTHORITY[“EPSG”,”32719″],
AXIS[“Easting”,EAST],
AXIS[“Northing”,NORTH]]
OGRFeature(estacion01_poly):0
OGRFeature(estacion01_poly):1
OGRFeature(estacion01_poly):2
OGRFeature(estacion01_poly):3

  • Con esta información de los EPSG ya somos capaces de reproyectar el tif desde “EPSG”,”4326″ a “EPSG”,”32719″ con el comando gdalwarp :

$ gdalwarp -s_srs EPSG:4326 -t_srs EPSG:32719 elqui-mosaico.tif elqui-mosaico_r.tif

  • Necesitamos generar el metadato de su Datum “EPSG”,”32719″, para nuestro shape de la divisoria de aguas de la cuenca estacion01_poly.shp con el comando epsg_tr.py :

$ epsg_tr.py -wkt 32719 > estacion01_poly.prj

  • Ya estamos listos para emascarar los datos que están fuera del polígono definido por nuestro shape en este caso el de a divisoria de aguas de la cuenca estacion01_poly.shp , con el comando gdalwarp -cutline :

$ gdalwarp -cutline estacion01_poly.shp elqui-mosaico_r.tif elqui-mosaico_re.tif

Con esto queda un archivo tif con las dimensiones del original, pero fuera de este la elevación está en cero.
Para recortar el archivo tif con las dimensiones del shape, podemos utilizar el comando gdalwarp -te si recordamos que las dimensiones del shape estacion01_poly.shp eran Extent: (392889.046537, 6679377.052909) – (415780.512166, 6711932.273227) , ¿No me creen? revisen más arriba lo que les dio ogrinfo:
$ gdalwarp -te 392889.046537 6679377.052909 415780.512166 6711932.273227 elqui-mosaico_re.tif elqui-mosaico_rer.tif

  • Para crear un shape llamado contour_250.shp con las curvas de nivel cada 250 mts. ocupamos el comando gdal_contour -a elev :

$ gdal_contour -a elev elqui-mosaico_rer.tif contour_250.shp -i 250.0

Luego visualizamos el resultado en Qgis:

Geotiff de la cuenca con curvas de nivel cada 250 mts.

BONUS: con el comando listgeo que nos entregan información interesante de los archivos Geotiff:
$ listgeo -no_norm ASTGTM_S29W068_dem.tif
$ listgeo -proj4 ASTGTM_S29W068_dem.tif
$ listgeo -t tabledir ASTGTM_S29W068_dem.tif

Referencias:
Manejando información SRTM con GDAL
Como instalar las herramientas acá utilizadas revisar: Software libre para Hidrología

Software libre para Hidrología


Actualmente me encuentro haciendo mi memoria de trabajo de título en “DESARROLLO DE UN MODELO DE PRONÓSTICO DE CAUDALES DE DESHIELO AFLUENTES AL EMBALSE PUCLARO UTILIZANDO HERRAMIENTAS DE SOFTWARE LIBRE”, para lo cual he tenido que instalar un entorno de desarrollo dentro de mi distribución favorita de GNU/Linux llamada Ubuntu.

Actualmente tengo instalada la versión 9.04 Jaunty Jackalope, pero si están comenzando les recomiendo que instalen la última versión desde ya, para comenzar a disfrutar lo más nuevo que ofrece el software libre.

Los paquetes que actualmente tengo instalados para trabajar son los siguientes:

Para cálculo numérico e interfaz gráfica:
– python (El lenguaje de desarrollo ágil y efectivo Python )
– ipython (Consola interactiva iPython)
– idle (Entorno de desarrollo integrado IDLE ), aunque Geany también me ha resultado cómodo de usar.
– python-scipy (Scientific tools for Python)
– python-numpy ( Numerical Python )
– python-rpy (Interfaz para trabajar con R, pero si instalan Ubuntu 9.10 recomiendo usar rpy2)
– python-matplotlib (Librería matplotlib para graficas 2D en python)
– python-tables (En caso que algún día tengamos muchos datos y desde muchas fuentes podremos probar PyTables Hierarchical database for Python based on HDF5)
– ( Tengo muchas ganas de instalar pytseries , pero sólo lo he visto disponible para Ubuntu 9.10 )

Para análisis estadístico:
– rkward (rkward Frontend de R)

Para trabajo con herramientas GIS:
qgis (Me ha sorprendido la versatilidad de este programa, nada que envidiar a los típicos ArcGis )
grass ( El clásico GIS opensource grass, lo bueno es activar los complementos de qgis para trabajar con los módulos de grass de una forma intuitiva y amigable)
marble (bajar imágenes atlas, el complemento ideal para obtener una imagen de la zona y georeferenciar algunos puntos )

Manipular GeoTiff:
–  gdal-bin (Permite además trabajar con elementos vectoriales)
–  geotiff-bin
–  python-imaging (Python Imaging Library PIL)
–  hdf4-tools (HDF command line utilities )
–  python-gdal Python (bindings to the Geospatial Data Abstraction Library)

Sigue leyendo