Obtención de atributos de un shapefile mediante el módulo ogr de Python

Siguiendo la tónica de artículos anteriores sobre el uso del módulo ogr de Python, en éste se van a obtener los atributos (campos de la tabla atributiva) de un shapefile conjuntamente con su geometría. Cooresponde a la página 41 del libro Python Geospatial Development. Como la parte inicial del código es coincidente con la de algunos precedentes se asigna a una función para poder ser utilizada más adelante; si es el caso. El código es el siguiente:

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

def despliegue():
	system('clear')
	path = '/home/zeito/Desktop/tl_2009_us_state'
	chdir(path)
	print "Archivos *.shp en %s" % path
	system('ls *.shp')
	nameshape = str(raw_input("Nombre del shapefile = ? "))
	shapefile = ogr.Open(nameshape)

	if shapefile is None:
		print "El shapefile no existe"
		exit(1)

	else:
		print "El shapefile se abrió satisfactoriamente"

	return shapefile

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

shapefile = despliegue()

que ejecutado en cónsola de bash produce la siguiente salida:

Archivos *.shp en /home/zeito/Desktop/tl_2009_us_state
tl_2009_us_state.shp
Nombre del shapefile = ? tl_2009_us_state.shp
El shapefile se abrió satisfactoriamente

e indica que el shapefile tl_2009_us_state.shp se abrió satisfactoriamente. Ahora, se van a obtener, sucesivamente, la capa (layer), los rasgos (features) y los atributos (items) del shapefile; sólo para el rasgo número 2 (Arizona). Estos últimos (atributos) se despliegan preliminarmente por pantalla para observar la estructura de lo que se imprime. El resto del código es el siguiente:

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

def despliegue():
	system('clear')
	path = '/home/zeito/Desktop/tl_2009_us_state'
	chdir(path)
	print "Archivos *.shp en %s" % path
	system('ls *.shp')
	nameshape = str(raw_input("Nombre del shapefile = ? "))
	shapefile = ogr.Open(nameshape)

	if shapefile is None:
		print "El shapefile no existe"
		exit(1)

	else:
		print "El shapefile se abrió satisfactoriamente"

	return shapefile

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

shapefile = despliegue()

layer = shapefile.GetLayer(0)

feature = layer.GetFeature(2)

atributos = feature.items()

print atributos

y su ejecución produce:

Archivos *.shp en /home/zeito/Desktop/tl_2009_us_state
tl_2009_us_state.shp
Nombre del shapefile = ? tl_2009_us_state.shp
El shapefile se abrió satisfactoriamente
{'DIVISION': '8', 'INTPTLAT': '+34.2099643', 'NAME': 'Arizona', 'STUSPS': 'AZ', 'FUNCSTAT': 'A', 'REGION': '4', 'LSAD': '00', 'AWATER': 1026257344.0, 'STATENS': '01779777', 'MTFCC': 'G4000', 'INTPTLON': '-111.6024010', 'STATEFP': '04', 'ALAND': 294207737677.0}

La salida de atributos (cuyo significado se encuentra en la Web del U.S. Census Bureau) es un diccionario pareado campos:valor y puede extraerse de manera más elegante mediante un loop de la siguiente manera:

.
.
.
namerasgo = feature.GetField("NAME")

print "Los campos del rasgo cuyo NAME es %s son:" % namerasgo

for campo,valor in atributos.items():

	print "%s = %s" % (campo,valor)
.
.
.

La geometría incluye adicionalmente esta porción de código:

.
.
.
geometria = feature.GetGeometryRef()

NombGeomet = geometria.GetGeometryName()

print "La geometria del rasgo %s es %s" % (namerasgo,NombGeomet)
.
.
.

cuya incorporación en el programa completo:

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

def despliegue():
	system('clear')
	path = '/home/zeito/Desktop/tl_2009_us_state'
	chdir(path)
	print "Archivos *.shp en %s" % path
	system('ls *.shp')
	nameshape = str(raw_input("Nombre del shapefile = ? "))
	shapefile = ogr.Open(nameshape)

	if shapefile is None:
		print "El shapefile no existe"
		exit(1)

	else:
		print "El shapefile se abrió satisfactoriamente"

	return shapefile

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

shapefile = despliegue()

layer = shapefile.GetLayer(0)

feature = layer.GetFeature(2)

atributos = feature.items()

namerasgo = feature.GetField("NAME")

print "Los campos del rasgo cuyo NAME es %s son:" % namerasgo

for campo,valor in atributos.items():

	print "%s = %s" % (campo,valor)

geometria = feature.GetGeometryRef()

NombGeomet = geometria.GetGeometryName()

print "La geometria del rasgo %s es %s" % (namerasgo,NombGeomet)

produce la siguiente salida:

Archivos *.shp en /home/zeito/Desktop/tl_2009_us_state
tl_2009_us_state.shp
Nombre del shapefile = ? tl_2009_us_state.shp
El shapefile se abrió satisfactoriamente
Los campos del rasgo cuyo NAME es Arizona son:
DIVISION = 8
INTPTLAT = +34.2099643
NAME = Arizona
STUSPS = AZ
FUNCSTAT = A
REGION = 4
LSAD = 00
AWATER = 1026257344.0
STATENS = 01779777
MTFCC = G4000
INTPTLON = -111.6024010
STATEFP = 04
ALAND = 294207737677.0
La geometria del rasgo Arizona es POLYGON

El significado de los campos se encuentra a continuación:

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