Encontrar los puntos más al norte, sur, este y oeste del feature de un shapefile mediante el módulo OGR de Python

Las páginas 43 y 44 del libro Python Geospatial Development presentan un código que, mediante el módulo OGR de Python, emplea una función recursiva que barre toda la geometría de un rasgo (feature) de un shapefile para obtener los puntos individuales que representan sus valores de mayor y menor latitud. No obstante, cabe señalar que en la función de la página 44 existe un error de imprenta que conllevó a una mala identación de los dos if que en ella se encuentra (son internos al for); lo cual se corrige en el código que se incluye más adelante.

Por otra parte, la función recursiva se amplía para incluir también en el diccionario la determinación de los puntos individuales que presentan la mayor y menor longitud del feature en el shapefile, lo que permite determinar dentro del cuerpo del programa las coordenadas de los vértices del cuadrilátero que abarca completamente el rasgo (que en este caso es el número 53, es decir, California).

El código propuesto es el siguiente:

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

from osgeo import ogr
from os import system,chdir

def findPoints(geometry, results):

	for i in range(geometry.GetPointCount()):
		x,y,z = geometry.GetPoint(i)
		
		if results['norte'] == None or results['norte'][1] < y:
			results['norte'] = (x,y)

		if results['sur'] == None or results['sur'][1] > y:
			results['sur'] = (x,y)

		if results['este'] == None or results['este'][0] < x:
			results['este'] = (x,y)

		if results['oeste'] == None or results['oeste'][0] > x:
			results['oeste'] = (x,y)

	for i in range(geometry.GetGeometryCount()):
		findPoints(geometry.GetGeometryRef(i), results)

system("clear")

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

chdir (path)

print "Shapefiles en %s" % path

system ("ls *.shp")

shapename = "tl_2009_us_state.shp"

shapefile = ogr.Open(shapename)

layer = shapefile.GetLayer(0)

numrasgo = int(raw_input("Número del rasgo = ? "))

feature = layer.GetFeature(numrasgo)

namefeature = feature.GetField("NAME")

print "Para el estado de %s" % namefeature

geometry = feature.GetGeometryRef()

results = {'norte' : None, 
		   'sur' : None,
		   'este' : None,
		   'oeste' : None}

findPoints(geometry, results)

print "El punto mas al norte es (%0.4f, %0.4f)" % results['norte']
print "El punto mas al sur es (%0.4f, %0.4f)" % results['sur']
print "El punto mas al este es (%0.4f, %0.4f)" % results['este']
print "El punto mas al oeste es (%0.4f, %0.4f)" % results['oeste']

n,N = results['norte']
s,S = results['sur']
O,o = results['oeste']
E,e = results['este']

print "Los vértices del cuadrilátero que abarca al rasgo %d son:" % numrasgo

print "Vértice 1: (%0.4f, %0.4f)\nVértice 2: (%0.4f, %0.4f)\nVértice 3: (%0.4f, %0.4f)\nVértice 4: (%0.4f, %0.4f)" % (O, N, E, N, E, S, O, S)

cuya salida, en cónsola de bash, es la siguiente:

Shapefiles en /home/zeito/Desktop/tl_2009_us_state/
tl_2009_us_state.shp
Número del rasgo = ? 53
Para el estado de California
El punto mas al norte es (-122.3782, 42.0095)
El punto mas al sur es (-117.2049, 32.5288)
El punto mas al este es (-114.1312, 34.2627)
El punto mas al oeste es (-124.4820, 40.4403)
Los vértices del cuadrilátero que abarca al rasgo 53 son:
Vértice 1: (-124.4820, 42.0095)
Vértice 2: (-114.1312, 42.0095)
Vértice 3: (-114.1312, 32.5288)
Vértice 4: (-124.4820, 32.5288)

Para verificar que el programa produce los resultados esperados los vértices, conjuntamente con el rasgo (individualizado del shapefile mediante Vectorial -> Herramientas de Gestión de datos -> Dividir Capa Vectorial en QGIS), se cargaron en una vista de QGIS. El shapefile de puntos fue pasado a polígono mediante el plugin Points2One y, tal como se observa en la siguiente imagen, circunscribe al estado de California.

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