Point shapefile, con librería ogr de Python, a partir de un archivo de texto plano

En un artículo anterior, se consideró un script de Python que crea un shapefile de puntos y los añade de manera interactiva a través de la cónsola. Sin embargo, si se dispone, por ejemplo, de gran cantidad de puntos es preferible que éstos sean añadidos a partir de un archivo de texto plano. El considerado (datos1) para probar el script fue el siguiente:

20076.783907 20311.382620
20252.801451 20433.252846
20419.119520 20463.071650
20459.532526 20628.198234
20447.667516 20650.203269
20543.655666 20623.512222
20528.621928 20558.220670
20633.378500 20551.075002
20623.429635 20171.705432
20848.667979 20004.109936
20916.476008 19695.501588
21006.620994 19222.932578
21117.249326 19116.260582
21157.424867 19042.070133
21427.706479 18619.992782
21711.369400 18134.324067
21140.450504 17875.180217
20915.216777 18197.905132
20814.265999 18325.588386
20504.089874 18730.825691
20205.874932 19118.796675
19885.441114 19535.423040
19998.989292 20001.727250

Como puede observarse, presenta por línea las coordenadas x,y separadas simplemente por un espacio. Fue obtenido de un artículo en el cual se convertían líneas de rumbo y distancia a coordenadas cartesianas. El script que se presenta a continuación usa inicialmente los módulos numpy y StringIO para convertir los datos a un formato de array. Posteriormente, el script tiene una lógica en la cual secuencialmente:

    1. Se establece el driver para el manejo de los shapefiles.
    2. Se crea una nueva fuente de datos y capa.
    3. Se añade un campo id al shapefile de salida.
    4. Dentro de un bucle se crea un nuevo objeto punto.
    5. Se crea un nuevo rasgo.
    6. Se añade el rasgo a la capa de salida.
    7. Se destruye la geometría y el rasgo.
    8. Finalmente, fuera del bucle, se cierra la fuente de datos.

El script propuesto es el siguiente:

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

from os import system,path

import numpy as np
from StringIO import StringIO 

from osgeo import ogr

system("clear")
 
name = str(raw_input("Nombre de shape a crear = ? "))

pfile=open(name,'r')
  
data=pfile.read()
  
pfile.close()

data=np.genfromtxt(StringIO(data)) #Se sobre entiende que los 
                                   #delimitadores son espacios
 
# Obtiene el driver
driver = ogr.GetDriverByName('ESRI Shapefile')

new_shape = name + '.shp'

# Crea una nueva fuente de datos y capa
if path.exists(new_shape):
    driver.DeleteDataSource(new_shape)
ds = driver.CreateDataSource(new_shape)

layer = ds.CreateLayer(new_shape, geom_type=ogr.wkbPoint)

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

for i in range(len(data)):

	# Crea un nuevo objeto punto
	point = ogr.Geometry(ogr.wkbPoint)
	point.AddPoint(data[i][0], data[i][1])

	# Obtiene el FeatureDefn para la capa de salida
	featureDefn = layer.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
	layer.CreateFeature(feature)
 
    # Destruye la geometría y el rasgo
	point.Destroy()
	feature.Destroy()

# Cierra la fuente de datos
ds.Destroy()

Su ejecución produce simplemente lo siguiente:

Nombre de shape a crear = ? datos1

y el resultado, datos1.shp, al visualizarlo en QGIS requiere que seleccionemos primeramente el sistema de coordenadas:

python1

cuya visualización en WGS 84/UTM 19N es la siguiente:

python2

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

6 respuestas a Point shapefile, con librería ogr de Python, a partir de un archivo de texto plano

  1. jerry dijo:

    Donde estan localizados los waypoints en el ejemplo. De acuerdo a las instrucciones es en Colombia?

  2. jerry dijo:

    Muchas gracias por contestar, soy de México y estoy estudiando QGis, su blog es de Venezuela?

    • Si. Por otra parte, QGIS es muy bueno si lo integras con GRASS y Sextante. Si finalmente logras también integrarlo con Python entonces no necesitas ArcGis. No obstante, la mejor integración se logra en ambientes tipo Unix como Linux porque puedes hacer scripts de bash tan poderosos como en Python sin necesidad de este último. Gracias por tu comentario.

      Saludos

      • jerry dijo:

        Gracias por contestar, estuve tratando de realizar reproyecciones de Conica de Lambert a UTM, y no pude con QGis 1.8, podria extenderse más en este tema. En cuanto a la utilización de GRASS en QGis, que texto me recomienda para entender mejor esta herramienta. Atentamente.

  3. Pingback: El módulo (librería) pyshp como soporte para leer/escribir archivos vectoriales en formato ESRI |

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