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

En un artículo precedente se consideró la manera de cómo extraer las coordenadas x, y de un shapefile tipo polígono mediante el módulo OGR de Python. Este sirve de base para transformarlo en un shapefile tipo punto adicionando los procedimientos de otro artículo en el cual se crea un Point shapefile vacío y se le añaden puntos. Además, se obtiene la referencia espacial del shapefile de origen (polígono) y se transfiere al shapefile de destino (puntos) simplemente exportándola como WKT a un archivo con extensión prj. Esto permite abrir directamente en QGIS los shapefiles creados sin que solicite la proyección respectiva (algo que no se había considerado en los artículos previos que usan el módulo OGR de Python).

El código a continuación se comenta para segregar los diferentes pasos de la transformación del shapefile tipo polígono cuenca.shp, almacenado en el directorio /home/zeito/Desktop/PRUEBAS_QGIS, para convertirlo en el shapefile de puntos cuenca_point.shp; ubicado en el mismo directorio anterior junto con el archivo cuenca_point.prj.

#!/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()

print spatialRef

#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_point.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_point.prj', 'w')
file.write(spatialRef.ExportToWkt())
file.close()

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

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

print "Coordenadas x,y"

for i in range(n_vertex):

	print x[i], y[i]

	# Crea un nuevo objeto punto
	point = ogr.Geometry(ogr.wkbPoint)
	point.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(point)
	feature.SetField('id', i+1)

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

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

La ejecución, en cónsola, del código anterior produce la siguiente salida:

Desde el directorio /home/zeito/Desktop/PRUEBAS_QGIS
El shapefile cuenca.shp se abrio satisfactoriamente
PROJCS["WGS_1984_UTM_Zone_19N",
    GEOGCS["GCS_WGS_1984",
        DATUM["WGS_1984",
            SPHEROID["WGS_1984",6378137,298.257223563]],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.017453292519943295]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-69],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]
numero de vertices =  15
El Shapefile cuenca_point.shp se creó exitosamente
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

Imprime más información que la necesaria pero es sólo con propósitos de verificación y puede ser eliminada si se desea. A continuación se tiene la imagen de la vista de QGIS que despliega los dos shapefiles. La tabla atributiva es la del shapefile creado (puntos).

Esta entrada fue publicada en Código Python, GDAL/OGR, SIG, Software Libre. Guarda el enlace permanente.

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

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

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