Economía para mejorar la calidad de vida

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

Filed under: chile | Tagged: gdal, gis, hydrology, ubuntu, ubuntu-cl |