Cómo extraer las coordenadas x, y de un polígono mediante el módulo OGR de Python

En un artículo precedente se señalaron las interesantes posibilidades de geoprocesamiento usando la librería OGR de Python. Sin embargo, se advertía el desconocimiento de los métodos, por ejemplo, para acceder directamente a las coordenadas de los vértices de un shapefile tipo polígono. Esto es relativamente sencillo para vectoriales tipo punto, tal como se trata en un artículo anterior, sin embargo, a medida que se profundiza en el curso de la Universidad de Utah se van descubriendo los procedimientos. Esto es importante porque ahorra el trabajo de transformar previamente el vectorial de polígono a punto para acceder a dicha información.

La clave del procedimiento para acceder a las coordenadas de los vértices está en comprender que el polígono está constituido por anillos (líneas) y la primera referencia geométrica que se solicita es con relación al polígono. Para acceder a las coordenadas de los vértices de los anillos hay que establecer una nueva referencia geométrica pero está vez con relación a los anillos. Por ejemplo, si se tiene un polígono con 5 anillos (referencias de 0 a 4) y se quieren las coordenadas de los vértices del tercer anillo la referencia debe apuntar al número 2.

En términos generales el procedimiento general involucra:

    1. Cambiarse al directorio de los shapefiles.
    2. Establecer el driver a usar.
    3. Abrir el archivo tipo polígono para solo lectura.
    4. Obtener la fuente de datos.
    5. Obtener el rasgo.
    6. Obtener la referencia geométrica para el poligono.
    7. Obtener la referencia geométrica para el anillo (linea).
    8. Se determina el numero de vertices de la linea.
    9. Cierra la fuente de datos y rasgos.

lo cual deriva en el siguiente código:

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

from os import chdir,system
from osgeo import ogr

system("clear")

path = '/home/zeito/Desktop/PRUEBAS_QGIS'

#Se cambia al directorio de los shapefiles
chdir(path)

#Establece el driver a usar
driver = ogr.GetDriverByName('ESRI Shapefile')

shapefile = 'cuenca.shp'

#Abre archivo tipo polígono para solo lectura
datasource = driver.Open(shapefile, 0)

print "Desde el directorio " + path

if datasource is not None:
	print "El shapefile " + shapefile + " se abrio satisfactoriamente"

else:
	print "El shapefile " + shapefile + " no existe"
	sys.exit(1)

#Obtiene la fuente de datos
layer = datasource.GetLayer()

#Obtiene el rasgo
feature = layer.GetFeature(0)

#Obtiene la referencia geometrica para el poligono
geometry = feature.GetGeometryRef()

#Obtiene la referencia geometrica para el anillo (linea)
vertice = geometry.GetGeometryRef(0)

#Se determina el numero de vertices de la linea
n_vertex = vertice.GetPointCount()

print "numero de vertices = ", n_vertex

print "Coordenadas x,y"

for i in range(n_vertex):
	print vertice.GetX(i), vertice.GetY(i)

# Cierra la fuente de datos y rasgos 
datasource.Destroy()
feature.Destroy()

La ejecución, en cónsola, para el vectorial cuenca.shp (polígono de un sólo anillo, id=1, referencia geométrica = 0 y 15 vértices) produce la siguiente salida:

Desde el directorio /home/zeito/Desktop/PRUEBAS_QGIS
El shapefile cuenca.shp se abrio satisfactoriamente
numero de vertices =  15
Coordenadas x,y
615534.042825 1051549.96936
626817.610648 1072537.40551
642388.934244 1087431.71504
668792.482951 1090591.11403
691359.618598 1089237.08589
710541.683898 1077050.83264
718214.510018 1065090.25075
716183.46781 1054258.02564
713926.754245 1041169.08696
709638.998472 1019053.29403
687974.548251 1008898.08299
644419.976453 1011380.46791
619144.784528 1021535.67895
612374.643834 1035527.30305
615534.042825 1051549.96936
Esta entrada fue publicada en Código Python, GDAL/OGR, SIG, Software Libre. Guarda el enlace permanente.

3 respuestas a Cómo extraer las coordenadas x, y de un polígono mediante el módulo OGR de Python

  1. Pingback: Transformar un shapefile tipo polígono en uno tipo punto mediante el módulo OGR de Python |

  2. pregunta: por ejemplo una parcela rectangular de 4 esquinas, y necesito las cuatro coordenadas que representan las esquinas.
    puedo con este script extraer los 4 puntos?

    • Si. El problema es que tus bindings de python, si es que usas Windows, apunten al sitio adecuado al emplear las librerías GDAL/OGR (el script está pensado para Linux donde el sistema resuelve las dependencias). Para ti es más fácil usar QGIS para resolver tu problema. Extraes los nodos de la parcela rectangular y luego, en ese vectorial de puntos, con la calculadora de campos, determinas las coordenadas. No es necesario script.

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