Transformar un shapefile tipo punto en uno tipo polígono mediante el módulo OGR de Python

En artículos anteriores se presentó código Python para, mediante el módulo OGR, transformar shapefiles de tipo polígono o línea en shapefiles tipo punto. En este se va a proceder en sentido inverso, es decir, transformar un shapefile tipo punto en otro tipo polígono. El procedimiento reutiliza gran parte del código ya presentado con anterioridad. Sin embargo, a pesar de que los rasgos en la capa (layer) pueden ser contados con la función GetFeatureCount(), se probó un método alternativo que implica el uso de un contador en un while para posteriormente acceder a las coordenadas en un bucle for con uso previo de ResetReading(); indispensable para una nueva lectura de rasgos (features).

Otro aspecto importante a ser considerado es que para obtener shapefiles de tipo polígono hay que crear primero anillos (líneas cerradas) que son incorporados posteriormente a la geometría del polígono. Habrá tantos anillos (grupos) como partes sencillas tenga el polígono.

El código propuesto es el siguiente:

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

from os.path import exists
from osgeo import ogr
from os import system
from os import chdir
import sys

#Borrar pantalla
system('clear') #En Win sería 'cls'

#Cambiarse al directorio donde se almacenan los shapefiles
path='/home/zeito/Desktop/PYTHON_GDAL/ospy_data1'
chdir(path)
print "Directorio donde se almacena el shapefile:\n", path

#Establecer el driver para manejar los vectoriales
driver = ogr.GetDriverByName('ESRI Shapefile')

#Abrir la fuente de datos (sólo lectura = 0)
shapefile = 'cuenca_point.shp'
datasource = driver.Open(shapefile, 0)

#Establecer excepción para apertura de archivo
if datasource is None:
	print "El Shapefile no existe"
	sys.exit(1)
else:
	print "El Shapefile " + shapefile + " se abrió exitosamente"

#Establecer la capa de la cual se extraeran los rasgos
layer = datasource.GetLayer()

#Establecer los campos que serán extraidos a partir de la capa
feature = layer.GetNextFeature()

#Establecer la referencia geometrica como rasgo
geometria=feature.GetGeometryRef()

#Obtiene la referencia espacial
spatialRef = layer.GetSpatialRef()

cont = 0

#Contando el numero de rasgos
while feature:
	cont += 1

	feature.Destroy()
	feature = layer.GetNextFeature()

layer.ResetReading()

feature = layer.GetNextFeature()

x = range(cont)
y = range(cont)

for i in range(cont):
	x[i] = geometria.GetX()
	y[i] = geometria.GetY()

	feature.Destroy()
	feature = layer.GetNextFeature()

datasource.Destroy()

# Crea una nueva fuente de datos y capa
fn = 'cuenca_pol.shp'

if exists(fn):
  driver.DeleteDataSource(fn)

outDS = driver.CreateDataSource(fn)

if outDS is None:
  print 'No puede crear el Shapefile' + fn
  sys.exit(1)

else:
  print 'El Shapefile ' + fn + ' se creó exitosamente'

# Copia la referencia espacial en el shapefile de salida
spatialRef.MorphToESRI()

file = open('cuenca_pol.prj', 'w')
file.write(spatialRef.ExportToWkt())
file.close()

# Para crear un poligono primero hay que crear un ANILLO

#Crea capa en el shapefile de salida tipo POLIGONO
outLayer = outDS.CreateLayer(fn, geom_type=ogr.wkbPolygon)

# Añade un campo id al shapefile de salida
fieldDefn = ogr.FieldDefn('id', ogr.OFTInteger)
outLayer.CreateField(fieldDefn)

# Crea un nuevo objeto tipo anillo
anillo = ogr.Geometry(ogr.wkbLinearRing)

# Añade los puntos al anillo
for i in range(cont):

	anillo.AddPoint(x[i], y[i])

#Cierra el anillo
anillo.CloseRings()

# Añade la geometría al poligono
polygon = ogr.Geometry(ogr.wkbPolygon)
polygon.AddGeometry(anillo)

# Obtiene el FeatureDefn para la capa de salida
featureDefn = outLayer.GetLayerDefn()

# Crea un nuevo rasgo
feature = ogr.Feature(featureDefn)
feature.SetGeometry(polygon)
feature.SetField('id', 1)

# Añade el rasgo a la capa de salida
outLayer.CreateFeature(feature)

# Destruye la geometría y el rasgo
anillo.Destroy()
polygon.Destroy()
feature.Destroy()

Ejecutado en cónsola, para el shapefile de puntos cuenca.shp, produce la salida siguiente en la obtención de cuenca_pol.shp:

Directorio donde se almacena el shapefile:
/home/zeito/Desktop/PYTHON_GDAL/ospy_data1
El Shapefile cuenca_point.shp se abrió exitosamente
El Shapefile cuenca_pol.shp se creó exitosamente

Para corroborar que el código funciona adecuadamente se tiene la siguiente imagen donde ambos shapefiles se despliegan en una vista de QGIS:

La tabla atributiva presentada en la imagen es la del shapefile tipo polígono.

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