Les comparto lo que ha sido en los últimos años mi acordeón (chuleta en España o cheat sheet en inglés) para GDAL/OGR.

Ojo: No están documentados ni todos los comandos ni todas las opciones, pero aquí están los que en algún momento he necesitado y me han sacado de apuros.

Última actualización Febrero 19 2011


Comandos para capas Raster (GDAL)

Obterner proyección, coordenadas extremas de la imagen, número de filas y columnas, tamaño del pixel, número de bandas y valores mínimos, máximos, promedios para cada banda de un archivo raster

gdalinfo -norat imagen_pancromatica.tif

Driver: GTiff/GeoTIFF
Files: imagen_pancromatica.tif
imagen_pancromatica.tif.aux.xml
Size is 11408, 7928
Coordinate System is:
PROJCS[“WGS 84 / UTM zone 13N”,
GEOGCS[“WGS 84”,
DATUM[“WGS_1984”,
SPHEROID[“WGS 84”,6378137,298.257223563,
AUTHORITY[“EPSG”,”7030″]],
AUTHORITY[“EPSG”,”6326″]],
PRIMEM[“Greenwich”,0],
UNIT[“degree”,0.0174532925199433],
AUTHORITY[“EPSG”,”4326″]],
PROJECTION[“Transverse_Mercator”],
PARAMETER[“latitude_of_origin”,0],
PARAMETER[“central_meridian”,-105],
PARAMETER[“scale_factor”,0.9996],
PARAMETER[“false_easting”,500000],
PARAMETER[“false_northing”,0],
UNIT[“metre”,1,
AUTHORITY[“EPSG”,”9001″]],
AUTHORITY[“EPSG”,”32613″]]
Origin = (631798.500000000000000,2294561.500000000000000)
Pixel Size = (3.000000000000000,-3.000000000000000)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  631798.500, 2294561.500) (103d44’2.56″W, 20d44’45.88″N)
Lower Left  (  631798.500, 2270777.500) (103d44’8.96″W, 20d31’52.38″N)
Upper Right (  666022.500, 2294561.500) (103d24’19.45″W, 20d44’36.03″N)
Lower Right (  666022.500, 2270777.500) (103d24’27.50″W, 20d31’42.65″N)
Center      (  648910.500, 2282669.500) (103d34’14.61″W, 20d38’14.52″N)
Band 1 Block=11408×1 Type=Byte, ColorInterp=Gray
Min=24.000 Max=255.000
Minimum=24.000, Maximum=255.000, Mean=89.825, StdDev=31.251
Metadata:
STATISTICS_MINIMUM=24
STATISTICS_MAXIMUM=255
STATISTICS_MEAN=89.824844767238
STATISTICS_MEDIAN=83
STATISTICS_MODE=75
STATISTICS_STDDEV=31.251114116524
STATISTICS_HISTONUMBINS=256
STATISTICS_HISTOMIN=0
STATISTICS_HISTOMAX=255
STATISTICS_SKIPFACTORX=9
STATISTICS_SKIPFACTORY=9
LAYER_TYPE=athematic
STATISTICS_HISTOBINVALUES=0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|81
|162|243|648|729|1215|2187|3645|4941|9234|11421|14823|20979|29808|44712|50139|67
149|87075|110646|131544|158922|197883|241947|282285|339390|391230|456516|507384|
590409|652536|702999|766017|847341|903312|953289|1015578|1087992|1171179|1219779
|1291383|1341522|1388583|1456218|1496880|1560870|1588248|1600479|1645434|1648836
|1673217|1677267|1719630|1686096|1698570|1670382|1669086|1615707|1604205|1574721
|1519398|1474929|1436049|1405512|1358856|1307745|1269594|1220265|1185192|1155141
|1136106|1074141|1044981|995409|997920|943893|923481|881928|872856|826119|812025
|770715|747630|720333|695061|679590|653427|633015|614871|579636|575019|551286|54
6750|517104|499608|485271|465669|442260|412857|399573|381510|379971|351621|34109
1|320193|309582|307233|292977|277344|260820|250290|233361|226800|219267|213921|2
02824|191808|186300|182736|175689|167508|160785|160299|151632|144990|134622|1367
28|133083|128304|119394|113238|111132|105624|111456|103761|96795|95337|95013|927
45|82377|83106|86346|79461|77598|78651|69741|66258|66906|64071|61398|63585|62694
|58482|59130|54189|53784|50382|49005|49653|46494|44226|42930|43902|41472|42525|3
6774|36045|39609|37584|37746|36531|29160|33453|28917|30618|28431|32157|28431|294
03|26325|28350|26244|24867|23814|24381|26973|24462|21141|22275|23166|20007|22680
|22923|22680|24948|24381|25272|28026|26487|25920|57996|93312|103194|59211|47871|
51597|39204|23247|22680|16200|13284|10530|6480|5022|5589|8829|4536|2997|2592|348
3|2187|2592|2511|2916|2511|2106|2349|2673|3402|2268|1863|2025|64881|

Pasar un archivo raster de un formato a otro

gdal_translate -of ECW imagen_1.img imagen_1.ecw

gdal_translate -of gtiff “C:\imagenes\imagen_1.img” “C:\proyecto\imagen_1.tif

Reproyectar un raster

gdalwarp -s_srs EPSG:4326 -t_srs EPSG:32613 ameca_4326.tif ameca32613.tif

Generar un mosaico

Generar un mosaico sencillo de unas pocas imágenes

gdal_merge -o mosaico.tifof GTiff imagen_1.img imagen_2.img imagen_3.img imagen_4.img

Generar un mosaico de muchas imágenes

Para esto se requiere un archivo .txt en donde se encuentren listadas las imágenes a mosaiquear

gdal_merge -o c:\proyecto\mosaico_grande.tif -of GTiff –optfile c:\imagenes\imagenes.txt

El archivo imagenes.txt tiene el nombre de la imágenes a mosaiquear incluyendo la extensión y la ruta si es necesario, una por renglon

c:\imagenes\imagen_1.img
c:\imagenes\imagen_2.img
c:\imagenes\imagen_3.img
c:\imagenes\imagen_4.img

Interpolar isolineas

gdal_contour -a elev mde.img curvas_nivel_10m.shp -i 10.0

Recortar una imagen con un polígono

Tenemos un raster que queremos cortar con una polígono, esto se consigue en dos pasos; en el primer paso generamos una  nueva imagen pero esta mantiene el extent de la imagen original

gdalwarp -cutline poligono.shp mosaico.tif recorte_paso1.tif

Esto nos da una imagen así:


En el segundo paso debemos recortar la imagen utilizando el extent del polígono, podemos obtener el extent del polígono utilizando

ogrinfo -so -al poligono.shp

INFO: Open of `poligono.shp’
using driver `ESRI Shapefile’ successful.

Layer name: poligono
Geometry: Polygon
Feature Count: 1
Extent: (-103.585645, 20.360664) – (-103.194984, 20.602870)
Layer SRS WKT:
(unknown)
ESTADO: String (255.0)
NOMBRE: String (255.0)
REGION: String (255.0)
Origen: String (255.0)
Clave: Real (32.0)
Area_km: Real (33.16)

De la información que nos arroja lo que utilizaremos será el extent

gdalwarp -te-103.585645, 20.360664-103.194984, 20.602870recorte_paso.tifrecorte_final.tif

Y el resultado final debe ser algo así:


Generar un sombreado (hillshade)

Trabajaremos con el MDE que recien cortamos, para hacer un MDE revise esta entrada.

gdaldem hillshade mde.tif -z 5 -s 111120 -az 315 sombreado.tif

Comandos para capas Vectoriales (OGR)

Pasar un archivo vectorial de un formato a otro

MapInfo (.tab) a Shape de ESRI (.shp)

ogr2ogr -f “ESRI Shapefile” salida.shp entrada.tab

Shape de ESRI (.shp) a Google Earth (.kml) *El .shp DEBE estar en EPSG:4326

ogr2ogr -f “KML” salida.kml entrada.shp

Shape de ESRI (.shp) a Google Earth (.kml) con cambio de proyección, en este ejemplo de UTM 13N WGS84 (EPSG:32613) a Geográficas WGS84 (EPSG:4326)

ogr2ogr -f “KML” -s_srs EPSG:32613 -t_srs EPSG:4326 salida.kml entrada.shp

Obtener geomatría, número de elementos, extent, proyección y descripción de la tabla de atributos del archivo vectorial

ogrinfo -so -al mapa.shp

Ejemplo:

ogrinfo -so -al mpios_jal_2009.shp

INFO: Open of `mpios_jal_2009.shp’
using driver `ESRI Shapefile’ successful.

Layer name: mpios_jal_2009
Geometry: Polygon
Feature Count: 125
Extent: (427469.290418, 2097415.178217) – (865450.642252, 2515726.071005)
Layer SRS WKT:
GEOGCS[“WGS84(DD)”,
DATUM[“WGS84”,
SPHEROID[“WGS84”,6378137.0,298.257223563]],
PRIMEM[“Greenwich”,0.0],
UNIT[“degree”,0.017453292519943295],
AXIS[“Geodetic longitude”,EAST],
AXIS[“Geodetic latitude”,NORTH]]
NOMBRE: String (255.0)
REGION: String (255.0)
Origen: String (255.0)
Clave: Real (32.0)
Area_km: Real (33.16)
fuente: String (255.0)

Consultar un mapa con SQL

ogrinfo mpios_jal_2009.shp -ro -geom=NO -sql
“select REGION,NOMBRE,Area_km from mpios_jal_2009
where REGION like ‘COSTA NORTE‘”

INFO: Open of `mpios_jal_2009.shp’
using driver `ESRI Shapefile’ successful.

Layer name: mpios_jal_2009
Geometry: Polygon
Feature Count: 3
Extent: (427469.290418, 2097415.178217) – (865450.642252, 2515726.071005)
Layer SRS WKT:
GEOGCS[“WGS84(DD)”,
DATUM[“WGS84”,
SPHEROID[“WGS84”,6378137.0,298.257223563]],
PRIMEM[“Greenwich”,0.0],
UNIT[“degree”,0.017453292519943295],
AXIS[“Geodetic longitude”,EAST],
AXIS[“Geodetic latitude”,NORTH]]
REGION: String (255.0)
NOMBRE: String (255.0)
Area_km: Real (33.16)
OGRFeature(mpios_jal_2009):98
REGION (String) = COSTA NORTE
NOMBRE (String) = TOMATL┴N
Area_km (Real) =             3150.7407004800002000

OGRFeature(mpios_jal_2009):99
REGION (String) = COSTA NORTE
NOMBRE (String) = CABO CORRIENTES
Area_km (Real) =             1465.4810887900001000

OGRFeature(mpios_jal_2009):104
REGION (String) = COSTA NORTE
NOMBRE (String) = PUERTO VALLARTA
Area_km (Real) =             1096.1276369200000000

Contar elementos de un mapa

ogrinfo mpios_jal_2009.shp-al -sql “select count (NOMBRE) from mpios_jal_2009

INFO: Open of `mpios_jal_2009.shp’
using driver `ESRI Shapefile’ successful.

Layer name: mpios_jal_2009
Geometry: Polygon
Feature Count: 1
Layer SRS WKT:
(unknown)
count_NOMBRE: Integer (0.0)
OGRFeature(mpios_jal_2009):0
count_NOMBRE (Integer) = 125

Estadísticas sencillas

Obtener la superficies mínima, maxima, promedio, la suma y la cuenta de los elementos en un mapa

ogrinfo mpios_jal_2009.shp-al -sql “select min(Area_km), max(Area_km), avg(Area_km), sum(Area_km), count(Area_km) from mpios_jal_2009

INFO: Open of `mpios_jal_2009.shp’
using driver `ESRI Shapefile’ successful.

Layer name: mpios_jal_2009
Geometry: Polygon
Feature Count: 1
Layer SRS WKT:
(unknown)
min_Area_km: Real (33.16)
max_Area_km: Real (33.16)
avg_Area_km: Real (33.16)
sum_Area_km: Real (33.16)
count_Area_km: Integer (0.0)
OGRFeature(mpios_jal_2009):0
min_Area_km (Real) =               81.9129196635999930
max_Area_km (Real) =             3150.7407004800002000
avg_Area_km (Real) =              633.8409782981970000
sum_Area_km (Real) =            79230.1222872746290000
count_Area_km (Integer) = 125

Extraer información de un archivo a otro nuevo

Extraer los elementos de un mapa mediante una consulta SQL a un mapa nuevo

ogr2ogr nuevo_mapa_X.shp -whereCAMPO = ‘X‘” mapa_completo.shp

Ejemplo:

ogr2ogr region_costa_norte.shp -where “REGION = ‘COSTA NORTE‘” mpios_jal_2009.shp

Extraer los elementos de un mapa a otro vía SQL aplicando un cambio de proyección, en este ejemplo de UTM 13N WGS84 (EPSG:32613) a Geográficas WGS84 (EPSG:4326)

ogr2ogr -s_srs EPSG:32613 -t_srs EPSG:4326 region_costa_norte_4326.shp -progress -whereREGION = ‘COSTA NORTE'” mpios_jal_2009.shp

El mismo caso pero aplicando un cambio de formato, en este ejemplo de Shape de ESRI (.shp) a Google Earth (.kml)

ogr2ogr -f “KML” -s_srs EPSG:32613 -t_srs EPSG:4326 region_costa_norte_4326.kml -progress -where “REGION = ‘COSTA NORTE'” mpios_jal_2009.shp

Extraer los elementos de un mapa que estén contenido dentro de otro mapa (clip); por ejemplo, extraer los puntos del mapa de localidades que estén dentro de la región Costa Norte

ogr2ogr locs_region_costa_norte.shp -progress -clipsrc region_costa_norte.shp locs_jal_2005.shp

Combinar (merge) archivos Shape de ESRI (.shp)

Primero hacemos una copia de cualquiera de los archivos .shp que deseamos combinar, en este caso tomare Costa Norte

ogr2ogr costa_completa.shp region_costa_norte.shp

ahora combinarémos region_costa_sur.shp con el mapa costa_completa.shp

ogr2ogr -update -append costa_completa.shp region_costa_sur.shp -nln costa_completa

 

Costa completa

Combinar muchas capas en una a través de un proceso por lotes (batch en windows)

Ejemplo en un catastro se tiene la información desagregada por manzanas y queremos generar mapas por zona catastral, los mapas de predios se nombran de acuerdo a la zona catastral donde están prd-01-020.shp son los predios de la manzana 020 de la zona catastral 01, entonces para generar un nuevo mapa que se llame predios01.shp hacemos lo siguiente:

for %q in (prd-01-*.shp) do ogr2ogr -update -append -f “ESRI shapefile” predios-01.shp %q -nln predios-01

Fuentes:

http://www.gdal.org/ogr_utilities.html

http://developmentseed.org/blog/2009/jul/30/using-open-source-tools-make-elevation-maps-afghanistan-and-pakistan

Manejando información SRTM con GDAL

http://linfiniti.com/2010/12/a-workflow-for-creating-beautiful-relief-shaded-dems-using-gdal/

http://www.gdal.org/ogr/drv_shapefile.html

Anuncios

2 comentarios sobre “Geo-Acordeón

  1. Creo que para recortar una imagen con un polígono es mucho más eficaz recortar primero el extent y luego hacer el recorte del polígono.

  2. Saludos ,

    Estoy convirtiendo un archivo DXF a un archivo SHP y de ahí a un archivo KML para luego superponerlo sobre google Earth. El caso es que consigo el paso de unos a otros gracias a la utilización de comandos de la librería GDAL-OGR.

    Todo sale perfecto hasta que superpongo el KML en google earth (el problema también ocurre al superponer el SHP en gvSIG sobre el castastro). El problema es que sale un poco desplazado al noreste respecto al mapa. He conseguido arreglar el problema mediante una herramienta (‘ProjectionUtility’ de ArcView) que realiza un geoproceso llamado: Geocentric[9603], ED_1950_TO_WGS_1984_13, POSC Code:8045.

    El problema es que necesito hacer esto de forma automática en código. No se si conoce la librería GDAL-OGR para darme alguna guía o si me recomienda alguna otra DLL.

    Utilizando el siguiente comando OGR (ogr2ogr -s_srs “+proj=utm +zone=30 +ellpse=intl +towgs84=-84,-107,-120,5,6,3,18 +units=m +no_defs” -t_srs “+proj=latlong +datum=WGS84” shapeOut.shp shapeIn.shp) consigo desplazar el SHP, pero ahora sale movido ligeramente hacia el suroeste.

    ¿Puede ayudarme? Muchas gracias.

    PD: conteste a mi correo si no le importa

Responder

Introduce tus datos o haz clic en un icono para iniciar sesión:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s