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

En el artículo anterior se estableció un procedimiento para transformar un shapefile tipo polígono en uno tipo punto y, en éste, el correspondiente a la transformación del mismo shapefile tipo polígono pero en uno tipo línea usando el módulo OGR de Python. El código incluye pocas variaciones y se presenta a continuación:

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

from os.path import exists
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 la referencia espacial
spatialRef = layer.GetSpatialRef()

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

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

#Obtiene la referencia 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

x = range(n_vertex)
y = range(n_vertex)

for i in range(n_vertex):
	x[i] = vertice.GetX(i)
	y[i] = vertice.GetY(i)

datasource.Destroy()
feature.Destroy()

# Crea una nueva fuente de datos y capa
fn = 'cuenca_line.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_line.prj', 'w')
file.write(spatialRef.ExportToWkt())
file.close()

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

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

# Crea un nuevo objeto linea
linea = ogr.Geometry(ogr.wkbLineString)

for i in range(n_vertex):

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

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

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

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

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

Su ejecución, en cónsola, produce la salida siguiente; donde se han excluido aquellos print que en el artículo anterior hacían referencia a las opciones denominadas con “propósito de verificación”.

Desde el directorio /home/zeito/Desktop/PRUEBAS_QGIS
El shapefile cuenca.shp se abrio satisfactoriamente
numero de vertices =  15
El Shapefile cuenca_line.shp se creó exitosamente

La imagen siguiente da cuenta de la visualización de los shapefiles involucrados lo que permite corroborar la utilidad del script:

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