Posteado por: rengifo | diciembre 17, 2010

Exploración y proyección de datos raster y de vectores usando GDAL o OGR


GDAL es una librería de abstracción de datos geoespaciales. OGR por su parte es una libreria para la manipulacion de abstraciones simples, como por ejemplo vectores. Ambas creadas por Frank Warmerdam

GDAL/OGR es soportada por varios tipos de software, entre los cuales se destacan, GRASS GIS , QGIS, MapServer, SAGA GIS, MapServer y MIRONE, entres otros. Por su flexibilidad GDAL/OGR soportan además una gran variedad de formatos, lo cual hace de GDAL/OGR  librerías con gran potencial bien sea en combinación con otros software o solas.

En esta oportunidad me gustaría introducirles a 3 funcionalidades, de hecho las más sencillas ya que la navidad se acerca:

  • Exploración de raster usando comandos con GDAL
  • Exploración de vector usando comandos con OGR
  • Reproyección de vectores y transformación de shapes a formato kml (google earth) usando comandos en (OGR).

Para empezar tienen que bajar la librería GDAL  para windows   o linux dependiendo del sistema operativo. Frank Warmerdam agrupó varios software incluyendo GDAL y OGR en una herramienta llamada  FWtools, que contiene entre otros OpenEV, MapServer, Proj4, OGDI y Python.

Para comenzar les recomiendo bajar FWtools  he instalarlo. Una vez hecho esto dos iconos van a aparecer en el escritorio (sistema operativo windows) (ver imagen 1)

Uno de ellos es el OpenEV  y el otro es la entrada ventana de  comandos (dos) para accesar directamente las herramientas de FWtools, es decir GDAL y OGR. Aqui supongo que la mayoría de ustedes ya se ha familiarizado con el uso del leguaje de programación ms-dos en windows. Si no, no hay problema, la sintásis que uso en este caso es  bastante sencilla.

Caso 1: Exploración de raster usando comandos con GDAL

Empecemos con la exploración de una imagen raster. En este caso, como casi siempre selecionaré el módelo (MED) de mi ciudad natal Mérida, Venezuela. Supongamos que recibes un MED y lo tratas de abrir en otro programa, por ejemplo SAGA GIS y cuando se abre hay algunos procesamientos que no pueden llevarse a cabo o que tardan demasiando tiempo y no sabes por que? (La ocurrencia de este caso se le agradece a uno de mis lectores)

Un punto de partida para saber que sucede es chequear el sistema de coordenadas o proyección de la imagen raster. La razón es que algunos algorítmos solo computan resultados en metros, pero no en minutos, grados y segundos.

Ahora, la pregunta  es como puedo chequear la estructura interna de una imagen raster si no tengo uno de estos software a la mano que cuestan tanto dinero? Bien, GDAL puede ser  de invaluable ayuda en esta situación.

Option 1:

Abrir OpenEV y visualizar la imágen o MED ( imágen 2) para ello  tienes file>open

Ya abierto puedes explorar las características de la imágen en la ventana layers, haces click con el botón derecho del ratón y se abre una ventana que muestra todas las propiedades de la imagen raster (imágen 3)

Option 2:

Abrir la ventana de comando (dos), navegar hacia donde tienes la imágen usando los comando cd( para navegar) y dir (para ver el contenido).

En mi caso mi directorio se ve asi:

Contenido C:\Documentos\rengifo\escritorio

28.07.2010  07:35             1 436 mirone.lnk

09.12.2010  09:27             2 066 Mozilla Firefox.lnk

20.02.2008  12:25               619 mstsc.exe.lnk

09.09.2008  14:09               441 MWSnap 3.lnk

30.11.2010  16:30           351 232 mypaper.mxd

28.05.2010  08:18               515 Nvu.lnk

17.12.2010  10:29             1 485 OpenEV_FW.lnk

11.08.2010  09:38               683 OpenJUMP.lnk

08.01.2009  16:18               524 OpenProj.lnk

11.06.2010  12:12         4 543 092 outgrass.asc

17.12.2010  11:30         5 554 335 outgrassev.asc

Una vez localizado el archivo en este outgrass.asc, escribo  el siguiente comando en la ventana de comando (dos), gdalinfo outgrass.asc. la respuesta que se verá como esta abajo me dice lo siguiente:

C:\Documents and Settings\RZO\Skrivebord>gdalinfo outgrass.asc

Driver: AAIGrid/Arc/Info ASCII Grid (TIPO DE FORMATO)

Files: outgrass.asc (TIPO DE ARCHIVO)

  • Size is 750, 563 (CANTIDAD DE CELDAS EN PLANO VERTICAL y HORIZONTAL, ES DECIR TOTAL)
  • Coordinate System is `’ (EL SISTEMA DE COORDENADAS NO ESTA PRESENTE)
  • Origin = (243120.000000000000000,973770.000000000000000)
  • Pixel Size = (60.000000000000000,-60.000000000000000) (TAMAÑO DE LAS CELDAS o PIXELS)
  • Corner Coordinates: (COORDENADAS DE LAS 4 ESQUINAS)
  • Upper Left  (  243120.000,  973770.000) (COORDENADAS SON EN METROS; LO QUE INDICA QUE POSIBLEMENTE QUE LA IMAGEN HA SIDO PROYECTADA A UTM)
  • Lower Left  (  243120.000,  939990.000)
  • Upper Right (  288120.000,  973770.000)
  • Lower Right (  288120.000,  939990.000)
  • Center      (  265620.000,  956880.000)
  • Band 1 Block=750×1 Type=Float32, ColorInterp=Undefined (LA IMAGEN SOLO TIENE UNA  BANDA; NO TIENE NINGUNA INTERPOLACION DE COLOR Y TIENE UNA MATRIX NUMÉRICA DE NUMEROS FLOTANTES DE DOBLE PRECISION, ES DECIR 32)
  • NoData Value=-9999  (LUGARES DONDE NO EXISTE DATOS EL NUMERO ES  DECIR  -9999)

Como se habran dado cuenta a la primera vista la línea de comandos pareciera ser mas difícil, en realidad es más rápida y directa. Ambas opciones generan los mismo resultados.

CASO 2: Exploración de vectores usando comandos con OGR

Dentro de  la misma linea  de comando puede navegar ahora hacia un archivo datos que contenga datos de contenga datos de  vectores en este caso en formato shape (.shp)

Usando el siguiente comando: ogrinfo -al -so sids.shp, puede obtener información de la misma manera como  explique anteriormente. Respecto a los parámetros –al y –so.  El primero indica que debe hacer un conteo (al) mientras que el segundo (so) le dice que lo debe resumir.

Digamos que tienes un vector y  quieres explorar su  estructura, en  este caso el archivo se llama sids.shp. Entonces escribo como lo explique en caso 1, dentro de la ventana de comandos lo siguiente, ogrinfo -al -so sids.shp

El resultado se vera como como esto

INFO: Open of `sids.shp’

using driver `ESRI Shapefile’ successful. (TIPO DE FORMATO)

Layer name: sids (NOMBRE DEL VECTOR)

Geometry: Polygon (GEOMETRIA)

Feature Count: 100 (CATIDAD DE ABSTRACCIONES)

Extent: (-84.323853, 33.881992) – (-75.456978, 36.589649) (LOCALIZACION EN LAT/LON)

Layer SRS WKT:

(unknown) (TIPO DESISTEMA  DE REFERENCIA DE COORDENADAS, EN ESTE CASO DESCONOCIDO)

AREA: Real (12.3) (AREA)

PERIMETER: Real (12.3) (PERIMETRO)

CNTY_: Real (11.0) (ATRIBUTOS)

CNTY_ID: Real (11.0)

NAME: String (32.0)

FIPS: String (5.0)

FIPSNO: Real (16.0)

CRESS_ID: Integer (3.0)

BIR74: Real (12.6)

SID74: Real (9.6)

NWBIR74: Real (11.6)

BIR79: Real (12.6)

SID79: Real (9.6)

Caso 3: Reproyección de vectores y transformación de shapes a formato kml (google earth) usando comandos en (OGR).

En este caso vamos a reproyectar un vector del sistema utm a lat/long de manera que lo podamos visualisar en google earth. Además, nos interesa que cuando le hagamos click en Google Earth,  nos muestre  algunas características  del vector, sin necesidad  de ir a las propiedades., por ejemplo.

Bueno el primer paso es seguir los pasos del caso 2 para explorar el archivo. Asi que es escribiendo el comando ogrinfo -al -so myshape.shp en la ventana de comandos  descubro la siguiente información:

C:\Documents and Settings\RZO\Skrivebord>ogrinfo -al -so myshape.shp

INFO: Open of `myshape.shp’

using driver `ESRI Shapefile’ successful.

Layer name: myshape

Geometry: Polygon

Feature Count: 6

Extent: (11502.760000, 6714159.481742) – (68935.690000, 6743239.280000)

Layer SRS WKT:

PROJCS["WGS_1984_UTM_Zone_33N",

GEOGCS["GCS_WGS_1984",

DATUM["WGS_1984",

SPHEROID["WGS_1984",6378137,298.257223563]],

PRIMEM["Greenwich",0],

UNIT["Degree",0.017453292519943295]],

PROJECTION["Transverse_Mercator"],

PARAMETER["latitude_of_origin",0],

PARAMETER["central_meridian",15],

PARAMETER["scale_factor",0.9996],

PARAMETER["false_easting",500000],

PARAMETER["false_northing",0],

UNIT["Meter",1]]

OBJECTID: Real (16.0)

Navn: String (10.0)

Fjordareal: Real (33.16)

Shape_area: Real (33.16)

Shape_len: Real (33.16)

Navn2: String (14.0)

FTEMA: Real (16.0)

VISNINGSTY: Real (16.0)

RESTRTYPE: Real (16.0)

KOMPONENT: String (13.0)

KOSTHOLDAR: String (59.0)

Sist vurde: Real (16.0)

OBJTYPE: String (18.0)

OPPHEVET: Real (16.0)

LINK: String (31.0)

Que el sistema de projección es WGS84 en la zona utm 33  (metros) y por las coordenadas todo  parece indicar que estamos en norte. Además puedo ver que el vector consiste de 6 partes, que es poligono sencillo y que tiene 14  columnas de atributos. Si embargo, aun no se de que se  trata el  tema. Asi que me voy a EV y abro el vector y me voy a tools>tabular shape grid y obtengo la table de los atributos (ver imágen 4)

Ahora se que se trata de un polígono que se encuentra en los fjordos en este caso Noruega, por los nombres y que mismo contiene información sobre la contaminacion de quimicos  dentro del fjordo. Tambien veo que tiene un link que me puede llevar a encontrar mas información.

Conociendo los atributos ahora ya se que atributos me gustaria mostrar cuando  alguien  visualize  el vector en google earth: Nombre tipo de quimico, ultimo analisis y fuente de mayor informacion. Como hacemos esto con OGR?

Bueno en primer lugar tengo que solucionar el problema de reproyección. Si no me será imposible  trasformar el shape a kml.

Para ello debo usar un comando que se llama ogr2ogr dentro de la ventana de comando y se hace de la siguiente manera: ogr2ogr -t_srs EPSG:4326 geomyshape.shp myshape.shp, donde myshape.shp es el archivo original y el resultante es el geomyshape.shp.

Si aplico el paso (comando) del caso 2 para explorar el vector, ogrinfo -al -so geomyshape.shp me dare cuenta que el geomyshape a cambiado el sistema de  coordenadas del utm a latitud y longitud (ver abajo).

C:\Documents and Settings\RZO\Skrivebord>ogrinfo -al -so geomyshape.shp

INFO: Open of `geomyshape.shp’

using driver `ESRI Shapefile’ successful.

Layer name: geomyshape

Geometry: Polygon

Feature Count: 6

Extent: (6.124072, 60.275328) – (7.141715, 60.586385)

Layer SRS WKT:

GEOGCS["GCS_WGS_1984",

DATUM["WGS_1984",

SPHEROID["WGS_1984",6378137,298.257223563]],

PRIMEM["Greenwich",0],

UNIT["Degree",0.017453292519943295]]

OBJECTID: Real (16.0)

Navn: String (10.0)

Fjordareal: Real (33.16)

Shape_area: Real (33.16)

Shape_len: Real (33.16)

Navn2: String (14.0)

FTEMA: Real (16.0)

VISNINGSTY: Real (16.0)

RESTRTYPE: Real (16.0)

KOMPONENT: String (13.0)

KOSTHOLDAR: String (59.0)

Sist vurde: Real (16.0)

OBJTYPE: String (18.0)

OPPHEVET: Real (16.0)

LINK: String (31.0)

Una vez completada esta parte nos podemos dedicar al segundo problema: Como transformar de shape a kml. Esto se  hace con el siguiente comando, ogr2ogr  -f “KML” -t_srs EPSG:4326 geomyshape.kml geomyshape.shp, donde geomyshape.kml es el nombre de producto final (ver imagen en google earth) Y  como lo habiamos pensado la zona se encuentra en el norte, es decir  en la costa oeste de Noruega (ver imágen 5)

Bueno estimados lectores de esta forma me despido y les deseo a todos una feliz navidad y próspero año Nuevo 2011!

About these ads

Responses

  1. Hola Rengifo, te cuento que me estoy mudando a linux en particular a Ubuntu 11.04 :-D . Estoy usando googleearth, GRASS GIS y hoy me anime a probar el FWTools y me parece fantastico este programa, estoy siguiendo este mini tutorial que hiciste y me resulto de mucha ayuda.
    Pdta.- Saga gis en linux es un quilombo instalar… el que viene por defecto es el 2.04 y la ultima version disponible es la 2.07, creo que seguire nomas con SAGA en windows pues en linux… no le entiendo nada de como instalar la ultima version jejeje. Un saludo!!!

    • Hola Diego como estás? .. Me alegra saber que diste el salto al Ubuntu. Al principio cuesta un poco adaptarse pero a medida que trabajas con este sistema operativo cada vez te dan menos ganas de usar windows. Respecto a SAGA GIS aun no he instalado la versión 2.0.7 ya que habia estado haciendo otras cosas. Pero en una mirada rápida.. parece que la mejor forma de instalarlo es usando PPA. Con lo cual puedes bajar los archivos ya compilados.
      Una vez usando el comando:
      sudo gedit /etc/apt/sources.list

      al ingresar tu clave se abre el la lista de recurso del ubuntu.. Alli le pegas los dos comandos que estan abajo.

      deb http://ppa.launchpad.net/johanvdw/saga-gis/ubuntu natty main
      deb-src http://ppa.launchpad.net/johanvdw/saga-gis/ubuntu natty main

      Anteriormente a ello tienes que haber bajado el PPA key..

      Finalmente le dice sudo get update para actualizar e instalar el software.
      Mas información hacer de los paquetes y como instalarlos la encuentras aca. Lamentablemente en inglés

      https://launchpad.net/~johanvdw/+archive/saga-gis
      suerte y saludos
      Rengifo

  2. Hola rengifo, una consulto como puedo reproyectar archivos DEM con FWTools. Yo descargo archivos de la nasa ASTER GDEM con resolucion a 30 metros, pero estos me vienen en grados, segundos… etc. quiero convertirlos a UTM como lo hago con este programita. Un saludo diego

    • Hola Diego,

      Me imagino que lo que tienes dos archivos por cada raster que has bajado: un ASTGTM_dem.tif y un ASTGTM_num.tif. El que nos interesa aqui es el ASTGTM_dem.tif. Como lo ves en el archivo el mismo es un tif, lo que hace la transformación bastante sencilla usando una utilidad quel GDAL que se GDALwarp, y que los tif contienen la cordenadas en el mismo archivo y por consiguientes no necesitas archivo extras como en el caso de los shapes.

      Si estas usando Linux y tienes el gdal instalado la sintáxis seria la siguiente:

      gdalwarp -s_srs “+proj=latlong +datum=WGS84″ -t_srs “+proj=utm +zone=48 +datum=WGS84″ -rb ge_20031018_last.tif ge_20031018_last_utm48.tif

      En el texto arriba:

      [-s_srs "+proj=latlong +datum=WGS84"] le dice al gdalwarp que la imagen inicial tienes este sistema de referencia espacial

      [-t_srs "+proj=utm +zone=48 +datum=WGS84"] le dice gdalwarp que reproyecte la imagen a UTM – Zone 48 – WGS sistema de coordenadas

      [-rb ] usa un algoritmo bilinear para la reproyección de la imágen. ( tienes otras opciones como por ejemplo bicubico [-rn] [-rc] [-rcs] es decir near (default), cubic, cubicspline, lanczos— generalmente el funciona bastante bien con el algoritmo bilinear

      [ge_20031018_last.tif] nombre del archivo/ imágen a transformar
      [ge_20031018_last_utm48.tif] imágen/archivo resultante despues de la transformación

      **** Si usas gdalwarp en windows (DOS) el nombre de imágen resultante viene primero y la imágen a transformar después. Si no vas a recibir un error donde dice que no encuentra la imagen original.

      Suerte.. y dejame saber si funcionó
      Saludos

  3. Hola rengifo de que tiempo. Una vez ingresada la siguiente entrada (por si acaso mi zona utm es la 20)

    gdalwarp -s_srs “+proj=latlong +datum=WGS84″ -t_srs “+proj=utm +zone=20 +datum=WGS84″ -rb dem.tif dem00.tif

    me sale el siguiente error en el gdalwarp:

    ERROR 1: Translating source or target SRS failed:
    “+proj=latlong

    Soy nuevo en gdal desde la linea de comandos y no entiendo el error, a que se refiere?

    Por otro lado, te comento que insatale SAGA-GIS 2.07 en ubuntu y por lo visto hay un error en los menues, el TEAM de saga deberia adaptar la ventana a la nueva ventana UNITY de ubuntu 11.04.

    Una consulta mas, me podrias indicar un programa para cortar un archivo raster?. Un saludo y segui adelante.

  4. Ya pille el error, es que copie directo desde tu post y me salieron las comillas en cursiva que no reconocia la consola… pero modificando me salio perfecto. MUCHAS GRACIAS RENGIFO!!! necesitaba aprender ese detallito para un trabajo de exposicion que estoy preparando, sera una ponencia sobre el uso de software libre en la ingenieria civil.

    Estare buscando informacion sobre como recortar un raster con gdal. Un saludo.

    • Me alegra que hayas encontrado el error… algunas veces cometer errores cuando tienes que usar comandos. Respecto a tu pregunta sobre como cortar un raster usando gdal/gdalwarp Desearía decir que hay un solución simple usando gdal. Pero para personas iniciandose en el uso de gdal creo es método es demasiado complicado. Ya que el mismo usa python-bindings. Acá un ejemplo
      http://linfiniti.com/2009/09/clipping-rasters-with-gdal-using-polygons/

      En general el script se veria algo asi:

      rtnY = geoMatrix[4]
      pixel = int((x – ulX) / xDist)
      line = int((ulY – y) / xDist)
      return (pixel, line)

      def histogram(a, bins=range(0,256)):
      “””
      Histogram function for multi-dimensional array.
      a = array
      bins = range of numbers to match
      “””
      fa = a.flat
      n = gdalnumeric.searchsorted(gdalnumeric.sort(fa), bins)
      n = gdalnumeric.concatenate([n, [len(fa)]])
      hist = n[1:]-n[:-1]
      return hist

      def stretch(a):
      “””
      Performs a histogram stretch on a gdalnumeric array image.
      “””
      hist = histogram(a)
      im = arrayToImage(a)
      lut = []
      for b in range(0, len(hist), 256):
      # step size
      step = reduce(operator.add, hist[b:b+256]) / 255
      # create equalization lookup table
      n = 0
      for i in range(256):
      lut.append(n / step)
      n = n + hist[i+b]
      im = im.point(lut)
      return imageToArray(im)

      # Load the source data as a gdalnumeric array
      srcArray = gdalnumeric.LoadFile(raster)

      # Also load as a gdal image to get geotransform (world file) info
      srcImage = gdal.Open(raster)
      geoTrans = srcImage.GetGeoTransform()

      # Create an OGR layer from a Field boundary shapefile
      field = ogr.Open(“%s.shp” % shp)
      lyr = field.GetLayer(shp)
      poly = lyr.GetNextFeature()

      # Convert the layer extent to image pixel coordinates
      minX, maxX, minY, maxY = lyr.GetExtent()
      ulX, ulY = world2Pixel(geoTrans, minX, maxY)
      lrX, lrY = world2Pixel(geoTrans, maxX, minY)

      # Calculate the pixel size of the new image
      pxWidth = int(lrX – ulX)
      pxHeight = int(lrY – ulY)

      clip = srcArray[:, ulY:lrY, ulX:lrX]

      # Create a new geomatrix for the image
      geoTrans = list(geoTrans)
      geoTrans[0] = minX
      geoTrans[3] = maxY

      # Map points to pixels for drawing the field boundary on a blank
      # 8-bit, black and white, mask image.
      points = []
      pixels = []
      geom = poly.GetGeometryRef()
      pts = geom.GetGeometryRef(0)
      for p in range(pts.GetPointCount()):
      points.append((pts.GetX(p), pts.GetY(p)))
      for p in points:
      pixels.append(world2Pixel(geoTrans, p[0], p[1]))
      rasterPoly = Image.new(“L”, (pxWidth, pxHeight), 1)
      rasterize = ImageDraw.Draw(rasterPoly)
      rasterize.polygon(pixels, 0)
      mask = imageToArray(rasterPoly)

      # Clip the image using the mask
      clip = gdalnumeric.choose(mask, (clip, 0)).astype(gdalnumeric.uint8)

      # This image has 3 bands so we stretch each one to make them
      # visually brighter
      for i in range(3):
      clip[i,:,:] = stretch(clip[i,:,:])

      # Save ndvi as tiff
      gdalnumeric.SaveArray(clip, “%s.tif” % output, format=”GTiff”, prototype=raster)

      # Save ndvi as an 8-bit jpeg for an easy, quick preview
      clip = clip.astype(gdalnumeric.uint8)
      gdalnumeric.SaveArray(clip, “%s.jpg” % output, format=”JPEG”

      Ahora si necesitas cortar una imagen de manera rectangular. Te recomiendo bajar el QGIS e installar un plugin que se llama Image cutter.

      Ahora bien si en realidad lo quieres hacer con GDAL entonces tendrias que instalar las siguientes librerias:

      Python Imaging Library
      Numpy

      De manera que puedas usar el siguiente script
      http://geospatialpython.com/2011/02/clip-raster-using-shapefile.html
      Saludos
      Rengifo

  5. estoy buscando información de como puedo convertir una capa de un formato con extensión .tab a .shp o a JSON. Un saludo

    • Hola Alejandro como estás? Me disculpas la tardanza. Mira creo que para lo que necesitas lo mejor es que uses el módulo ogr2ogr de la librería GDAL. La sintáxis sería algo como así:

      ogr2ogr -f geoJSON copia.json original.shp

      Después de haber bajado el GDAL navegas a donde se encuentran los datos en linea de comando (en linux el shell y en windows el CMD:

      *llamas el modulo ogr2ogr
      *defines el formato al que quieres transformar -f geoJSON
      *defines el nombre del la copia a crear en geojson por ejemplo xxx.geojson
      * escribres el nombre del shp a transformar
      Solo lo he probado con shapes pero deberia funcinar con tab ya que GDAL tambien soporta el formato de mapInfo
      Saludos

    • Tambien hay una posibilidad que igualmente usa gdal en el fondo. Puedes installar gdal y despues bajar el QGIS. Con el mismo puedes transformar entre diferente formatos..!


Deja un comentario

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

Categorías

Seguir

Recibe cada nueva publicación en tu buzón de correo electrónico.

%d personas les gusta esto: