Análisis preliminar de shapefiles usando módulo OGR de Python

En un artículo precedente se inició el estudio del libro “Python Geospatial Development” de Erik Westra mostrando un ejemplo en el cual se determinaba la distancia entre dos puntos de la superficie terrestre cuyas coordenadas estaban en geográficas (grados decimales) y considerando el elipsoide WGS84 como referencia. En este artículo se pretende desarrollar de manera comprensiva el código de la página 40 al modo de Learn Python The Hard Way. La razón es para evitar un simple “copy and paste” del código para “observar” que funciona y poder escribirlo razonadamente sin uso de material de apoyo.

El shapefile usado para el ejercicio fue bajado de aquí:

http://www2.census.gov/geo/tiger/TIGER2009/tl_2009_us_state.zip

y lo descomprimí en el escritorio de mi home de usuario (Linux Debian). Por tanto, la ruta al shapefile, en mi caso, es:

/home/zeito/Desktop/tl_2009_us_state/tl_2009_us_state.shp

Como algunas aspectos del código parecen “surgidos de la manga”, es obvio que para escribirlo de esa manera se tenía un conocimiento previo de que contenía. Para ello, el shapefile se desplegó en una vista de QGIS conjuntamente con su tabla atributiva y se expone en la imagen siguiente:

Corresponde a un shapefile tipo polígono que representa a Estados Unidos (USA) con una tabla atributiva que contiene un campo denominado NAME y que refiere el nombre de sus estados, estados libres asociados, etc.

El código de la página 40, a grandes rasgos, pretende establecer cuántos layers tiene el shapefile y cuál es su referencia espacial, cuantos features tiene cada layer y desplegar el nombre del campo NAME para cada feature.

En mi caso, yo comienzo mis scripts en Linux Debian con las siguientes dos líneas:

#!/usr/bin/ebv python
# -*- coding: utf-8

La segunda me permite escribir en los comentarios y en las salidas caracteres especiales como acentos y eñes; por ejemplo. Otra cosa que me agrada es limpiar la pantalla y como tengo todos los scripts de python ubicados en un mismo sitio, requiero entonces mudarme al directorio en el cual se encuentra el shape. Una vez allí sería útil, como verificación, hacer un listado de los shapefiles del directorio. Todo esto se logra con comandos (system, chdir) propios del módulo os. Esto daría lugar al siguiente código:

#!/usr/bin/ebv python
# -*- coding: utf-8

from os import system,chdir

system("clear")

path = "/home/zeito/Desktop/tl_2009_us_state/"

chdir(path)

print "Los shapefiles en %s son:" % path

system("ls *.shp")

cuya ejecución, en cónsola, produce lo siguiente:

Los shapefiles en /home/zeito/Desktop/tl_2009_us_state/ son:
tl_2009_us_state.shp
zeito@debian:~$ 

Para poder acceder a la información del shapefile es necesario abrirlo con el método Open del módulo ogr de python (englobado en python-gdal). Para instalarlo en Debian, como superusuario:

aptitude install python-gdal

Para invocarlo en el script e incluir opciones de verificación de existencia del shapefile tendríamos:

#!/usr/bin/ebv python
# -*- coding: utf-8

from os import system,chdir
import sys
from osgeo import ogr

system("clear")

path = "/home/zeito/Desktop/tl_2009_us_state/"

chdir(path)

print "Los shapefiles en %s son:" % path

system("ls *.shp")

NameShape = str(raw_input("Nombre del shapefile = ? "))

shapefile = ogr.Open(NameShape)

if shapefile is None:
	print "El shapefile %s no existe" % NameShape
	sys.exit(1)

else:
	print "El shapefile %s se abrió satisfactoriamente" % NameShape

Ahora, una vez abierto el shapefile, se tiene que contar sus layers con GetLayerCount y obtener cada una de ellas dentro de un bucle, su referencia espacial, así como, contar su número de rasgos con GetFeatureCount. Esto añadiría la siguiente porción de código:

.
.
.
for NumCapas in range(numcapas):

	layer = shapefile.GetLayer(NumCapas)
	numrasgos = layer.GetFeatureCount()
	RefSpatial = layer.GetSpatialRef().ExportToProj4()	

	print "La referencia espacial de la capa %d es %s" % (NumCapas,RefSpatial)
	print "El número de rasgos de la capa %d es %d" % (NumCapas, numrasgos)
.
.
.

Finalmente, se tendría que añadir dentro de un bucle que barre todos los rasgos, la rutina para obtenerlos a partir de la capa (GetFeature) y a partir de los rasgos (features) obtener el campo “NAME” (GetField) e imprimirlo. El código quedaría finalmente así:

#!/usr/bin/ebv python
# -*- coding: utf-8

from os import system,chdir
import sys
from osgeo import ogr

system("clear")

path = "/home/zeito/Desktop/tl_2009_us_state/"

chdir(path)

print "Los shapefiles en %s son:" % path

system("ls *.shp")

NameShape = str(raw_input("Nombre del shapefile = ? "))

shapefile = ogr.Open(NameShape)

if shapefile is None:
	print "El shapefile %s no existe" % NameShape
	sys.exit(1)

else:
	print "El shapefile %s se abrió satisfactoriamente" % NameShape

numcapas = shapefile.GetLayerCount()

for NumCapas in range(numcapas):

	layer = shapefile.GetLayer(NumCapas)
	numrasgos = layer.GetFeatureCount()
	RefSpatial = layer.GetSpatialRef().ExportToProj4()	

	print "La referencia espacial de la capa %d es %s" % (NumCapas,RefSpatial)
	print "El número de rasgos de la capa %d es %d" % (NumCapas, numrasgos)

for NumRasgo in range(numrasgos):

	feature = layer.GetFeature(NumRasgo)

	NombreCampo = feature.GetField("Name")

	print "El nombre del campo NAME en el rasgo %d es %s" % (NumRasgo, NombreCampo)

que ejecutado en cónsola de bash produciría la siguiente salida:

Los shapefiles en /home/zeito/Desktop/tl_2009_us_state/ son:
tl_2009_us_state.shp
Nombre del shapefile = ? tl_2009_us_state.shp
El shapefile tl_2009_us_state.shp se abrió satisfactoriamente
La referencia espacial de la capa 0 es +proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs 
El número de rasgos de la capa 0 es 56
El nombre del campo NAME en el rasgo 0 es American Samoa
El nombre del campo NAME en el rasgo 1 es Nevada
El nombre del campo NAME en el rasgo 2 es Arizona
El nombre del campo NAME en el rasgo 3 es Wisconsin
.
.
.
El nombre del campo NAME en el rasgo 51 es Illinois
El nombre del campo NAME en el rasgo 52 es New Jersey
El nombre del campo NAME en el rasgo 53 es California
El nombre del campo NAME en el rasgo 54 es Ohio
El nombre del campo NAME en el rasgo 55 es Texas
Esta entrada fue publicada en Código Python, GDAL/OGR, SIG, Software Libre. Guarda el enlace permanente.

Responder

Por favor, inicia sesión con uno de estos métodos para publicar tu comentario:

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