Copia parcial de un vectorial de puntos con generación de nuevos campos, incluyendo coordenadas, con librería OGR de python

En un artículo precedente, se consideró la impresión de campos en archivo de texto, incluyendo coordenadas, con la librería OGR de Python. En este se va a señalar cómo copiar parcialmente la información de un vectorial de puntos existente a otro que se ha creado para tal fin pero en el cual también se han añadido nuevos campos que no existían en su predecesor.

Para ello, con el shapefile de origen, los procedimientos siguen una rutina basada más o menos en estos pasos:

1. Establecer el driver para manejar los vectoriales (‘ESRI Shapefile’).
2. Abrir la fuente de datos (‘sites.shp’).
3. Establecer excepción para apertura de archivo.
4. Establecer la capa de la cual se extraerán los rasgos.

sin embargo, para el nuevo shapefile los procedimientos requieren la creación de un nuevo layer y la copia o la creación de los nuevos rasgos, campos y la geometría. Como clave en la creación de los campos está el índice FieldDefn. Para la copia, el FieldDefn se define a través de los rasgos (features) usando la función miembro GetFieldDefnRef. Para la creación de nuevos campos la definición pasa entonces por la función ogr.FieldDefn() cuyo argumento depende del tipo de dato (entero, real, string). El script propuesto es el siguiente:

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

# Importando módulos
import os, sys
from osgeo import ogr
from math import sqrt

os.system("clear")

# Estableciendo el directorio de trabajo
os.chdir('/home/zeito/Desktop/PYTHON_GDAL/ospy_data1')

# Obteniendo el driver para manejar shapefiles
driver = ogr.GetDriverByName('ESRI Shapefile')

# Abriendo la fuente de datos (modo sólo lectura) y abriendo la capa
shape_input = 'sites.shp'
inDS = driver.Open(shape_input, 0)

if inDS is None:
  print 'El Shapefile ' + shape_input + ' no existe'
  sys.exit(1)

else:
	print 'El Shapefile ' + shape_input + ' se abrió satisfactoriamente'

inLayer = inDS.GetLayer()

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

if os.path.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'

outLayer = outDS.CreateLayer(fn, geom_type=ogr.wkbPoint)

# Obtiene los FieldDefn's para los campos id y cover desde el shapefile de entrada
feature = inLayer.GetFeature(0)
idFieldDefn = feature.GetFieldDefnRef('id')
coverFieldDefn = feature.GetFieldDefnRef('cover')

# Crea los nuevos campos id y cover en el shapefile de salida
outLayer.CreateField(idFieldDefn)
outLayer.CreateField(coverFieldDefn)

#Crear el FieldDefn para el campo x,y, distancia en el shapefile de salida
id_oldFieldDefn = ogr.FieldDefn('id_old', ogr.OFTInteger)
xFieldDefn = ogr.FieldDefn('x', ogr.OFTReal)
yFieldDefn = ogr.FieldDefn('y', ogr.OFTReal)
distFieldDefn = ogr.FieldDefn('dist', ogr.OFTReal)

# Crear nuevos campos x, y, id_old y distancia en el shapefile de salida
outLayer.CreateField(id_oldFieldDefn)
outLayer.CreateField(xFieldDefn)
outLayer.CreateField(yFieldDefn)
outLayer.CreateField(distFieldDefn)

if idFieldDefn is not None:
	print "El campo id en " + fn + " se creó exitosamente"

if coverFieldDefn is not None:
	print "El campo cover en " + fn + " se creó exitosamente"

if id_oldFieldDefn is not None:
	print "El campo id_old en " + fn + " se creó exitosamente"

if xFieldDefn is not None:
	print "El campo x en " + fn + " se creó exitosamente"

if yFieldDefn is not None:
	print "El campo y en " + fn + " se creó exitosamente"

if distFieldDefn is not None:
	print "El campo dist en " + fn + " se creó exitosamente"

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

# bucle a traves de los rasgos de entrada
inFeature = inLayer.GetNextFeature()

cont = 1

while inFeature:

	# Obtiene el atributo cover para el rasgo de entrada
	cover = inFeature.GetField('cover')

	# Verifica si cover == trees
	if cover == 'trees':

		# Crea un nuevo rasgo
		outFeature = ogr.Feature(featureDefn)

		# Establece la geometria
		geom = inFeature.GetGeometryRef()
		outFeature.SetGeometry(geom)
		x = geom.GetX()
		y = geom.GetY()

		if cont == 1:

			x_1 = x
			y_1 = y

		d = sqrt((x-x_1)*(x-x_1) + (y-y_1)*(y-y_1))

		# Establece los atributos
		id = inFeature.GetField('id')
		outFeature.SetField('id', cont)
		outFeature.SetField('id_old', id) 
		outFeature.SetField('cover', cover)
		outFeature.SetField('x', x)
		outFeature.SetField('y', y)
		outFeature.SetField('dist', d)

		cont += 1
 
		# Añade el rasgo a la capa de salida
		outLayer.CreateFeature(outFeature)

		# Destruye el rasgo de salida
		outFeature.Destroy()


	# Destruye el rasgo de entrada y se obtiene uno nuevo
	inFeature.Destroy()
	inFeature = inLayer.GetNextFeature()

# Cierra las fuentes de datos
inDS.Destroy()
outDS.Destroy()

que se aplicó al shapile sites.shp, del curso ya referido en el artículo anterior, para seleccionar y copiar del campo cover en el nuevo shapefile (arboles_point.shp) los puntos correspondientes a tree (árboles). La ejecución del script produce, en cónsola, este resultado:

El Shapefile sites.shp se abrió satisfactoriamente
El Shapefile arboles_point.shp se creó exitosamente
El campo id en arboles_point.shp se creó exitosamente
El campo cover en arboles_point.shp se creó exitosamente
El campo id_old en arboles_point.shp se creó exitosamente
El campo x en arboles_point.shp se creó exitosamente
El campo y en arboles_point.shp se creó exitosamente
El campo dist en arboles_point.shp se creó exitosamente

que da cuenta de los campos que fueron copiados y/o creados. En la imagen siguiente se tiene el shapefile creado con 11 puntos de los 41 totales en sites.shp y mostrando la correspondiente tabla atributiva.

Este shapefile está proyectado para la zona 12 UTM y abarca los estados Utah (principalmente) e Idaho; tal como se ve en Google Earth:

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