Hace unas semanas un amigo del proyecto, Gustavo Agüero , publicó un post
en su blog en el que explicaba cómo generar una poligonal en gvSIG a partir
de un rumbo o derrotero. Me pareció un ejercicio muy interesante y con su
colaboración implementamos un script que realizase los mismos pasos en gvSIG
2.0. El resultado es el que explicamos en esta entrada.

Quiero aclarar que todo el mérito de este ejerccio es de Gustavo y los errores
que hay son míos por no seguir correctamente sus indicaciones y por mi
ignorancia. Todo lo que veáis comentadlo y lo corregiré 😉

Partimos de un rumbo que presenta los dotos con la estructura LINEA, AZIMUT,
DISTANCIA, siendo AZIMUT el ángulo de una dirección contado en el sentido de
las agujas del reloj a partir del norte geográfico (descargar el rumbo.jpg).

Lo primero que debemos hacer es pasar este archivo a un formato que luego
seamos capaces de leer desde un script, este archivo que hemos creado es un
archivo csv al que hemos llamado rumbo.csv (descargar el rumbo.csv).

Una vez que tenemos nuestro archivo pensemos qué necesitamos para obtener
su contenido en un shape file. Obviamente necesitaremos un shape file para
guardar el resultado por lo que también necesitamos crear las entidades que
representen nuestros datos para
incluirlas el el shape. Para crear estas entidades necesitaremos transformar
las coordenadas polares (Azimut, distancia) en rectangulares (X,Y).

Un par de comentarios más, Gustavo en su post recomienda hacer una
representación de los datos. A mí me resultó de mucha ayuda.

En el script vamos a crear un shape de tipo LINEAS y otro de tipo
POLIGONOS. En el de tipo LINEAS, guardaremos los datos del rumbo, los datos
de la transformació de las coordenadas, y las coordenadas resultantes. En
el de POLIGONOS guardaremos un identificador de la entidad y un campo de
tipo texto.

Partimos de un gvSIG desktop 2.0 (en mi caso el build que he usado es
el 2060 RC1) con la última versión de la extensión de scripting instalada
(actualmente es la 36).

Una vez que tenemos abierto el editor de scripts (barra de menús
Herramientas/Scripting/Scripting Composer) creamos un script nuevo y empezamos
a escribir nuestro script.

Lo primero que vamos a hacer es crear las capas de salida usando la función
createShape del módulo gvsig. La sintaxis de esta función nos dice que
necesita una definición de datos para las entidades, una ruta que es el lugar
dónde debe almacenarse el shape file que creemos, la proyección del shape file,
y el tipo de geometría que va a contener. La definición de los datos la
creamos usando la función createSchema del mismo módulo gvsig .

El código con comentarios sería:

import gvsig
import geom

def main():

'''
Crea una capa de poligonos con el siguiente modelo de datos:
    - 'ID', Integer
    - 'OBSERVACIONES', string
    - 'GEOMETRY', Geometry

Crea una capa de lineas con el siguiente modelo de datos:
    - "ID", Integer
    - "LINEA", string
    - "GRADOS", long
    - "MINUTOS",long
    - "DISTANCIA", double
    - "RADIAN", double
    - "X", double
    - "Y", double
    - "GEOMETRY", Geometry            
'''

#Establecemos la proyeccion, es la misma para las 2 capas
CRS="EPSG:32617"

#Crea el objeto que representa el modelo de datos del shape de poligonos
schema_poligonos = gvsig.createSchema() 

#Inserta los campos del modelo de datos
schema_poligonos.append('ID','INTEGER', size=7, default=0) 
schema_poligonos.append('OBSERVACIONES','STRING', size=200, default='Sin modificar') 
schema_poligonos.append('GEOMETRY', 'GEOMETRY')

#Establecemos la ruta. Recuerda cambiarla!!
ruta='/tmp/rumbo-poligonos.shp'

#Creamos el shape
shape_poligonos = gvsig.createShape(
        schema_poligonos, 
        ruta,
        CRS=CRS,
        geometryType=geom.SURFACE
    )

#Creamos la capa de lineas

#Crea el objeto que representa el modelo de datos del shape de lineas
schema_lineas = gvsig.createSchema() 

#Inserta los campos del modelo de datos
schema_lineas.append('ID','INTEGER', size=7, default=0) 
schema_lineas.append('LINEA','STRING',size=50,default='')
schema_lineas.append('GRADOS','LONG', size=7, default=0) 
schema_lineas.append('MINUTOS','LONG', size=7, default=0)
schema_lineas.append('SEGUNDOS','LONG', size=7, default=0) 
schema_lineas.append('DISTANCIA','DOUBLE', size=20, default=0.0, precision=6) 
schema_lineas.append('RADIAN','DOUBLE', size=20, default=0.0, precision=6) 
schema_lineas.append('X','DOUBLE', size=20, default=0.0, precision=6) 
schema_lineas.append('Y','DOUBLE', size=20, default=0.0, precision=6) 
schema_lineas.append('GEOMETRY', 'GEOMETRY')

#Establecemos la ruta. Recuerda cambiarla!!
ruta='/tmp/rumbo-lineas.shp'

#Creamos el shape
shape_line = gvsig.createShape(
    schema_lineas, 
    ruta,
    CRS=CRS,
    geometryType=geom.MULTILINE
)

Bien, ya tenemos nuestras capas de salida, veamos cómo obtener los datos del
archivo csv que hemos creado. Python tiene un módulo que nos permite
manejar archivos csv de forma muy cómoda ( python csv ). El código para
leer el archivo csv sería:

import csv
import os.path

#Ruta donde esta el archivo csv. Recuerda cambiarlo!!
csv_file = '/tmp/rumbo.csv'

#Comprueba que el archivo exista en la ruta
if not os.path.exists(csv_file):
  print "Error, el archivo no existe"
  return

#Abrimos el archivo en modo lectura
input = open(csv_file,'r')

#Creamos un objeto reader del modulo csv
reader = csv.reader(input)

#Leemos el archivo
for row in reader:
    #imprimimos la linea
    print ', '.join(row)

Ya sabemos crear shapes y leer archivos csv , el siguiente paso sería crear
las entidades, para ello debemos transformar las coordenadas polares a
rectangulares y es en este punto donde para mí las cosas se vuelven un tanto
oscuras por lo que perdonad los errores, que los habrá, que podáis encontrar.
Para realizar la transformación convertiremos el ángulo decimal a un ángulo en
radian. Gustavo en su post cómo generar una poligonal entre los archivos
para descargar nos facilita un pdf con una explicación de los cálculos que ha
realizado, si alguno tiene dudas recomiendo que lo lea, aclarar únicamente que
usaremos el módulo de python math para usar las funciones del seno y del
coseno, necesarias para realizar la transformación. Además, nuestro
rumbo sabemos que utiliza coordenadas relativas, por lo que vamos a necesitar
un punto de partida y el código para calcular las nuevas coordenadas a partir
de las anteriores. Nuestro punto de partida es (361820.959424, 1107908.627000).
El codigo podría ser:

"""
Convert decimal degrees, minutes, seconds to radian
"""
import math
#Definimos el valor de la constante 'pi' que vamos a usar
PI = 3.1415926

x_ant = 361820.959424
y_ant = 1107908.627000

#Aplicamos la formula y obtenemos el angulo en radianes
angulo = (degrees+(minutes/60.0)+(seconds/3600.0))*PI/180.0

#Convertimos las coordenadas polares en rectangulares

x_new = radio * math.sin(angulo)
y_new = radio * math.cos(angulo)

#Al ser coordenadas relativas sumanos los valores de las coordenadas anteriores
x = x_new + x_ant
y = y_new + y_ant

La última pieza de nuestro puzzle es cómo creamos las entidades y sus
geometrías. Las capas que hemos creado anteriormente shape_line y
shape_poligonos tienen un método llamado append que nos permite añadir los
fenómenos o entidades a la capa. Este método acepta un diccionario,
python dict que usa como claves los campos que hemos definido en la
estructura de datos de la capa cuando la hemos creado, y los valores son los
que debe tener la entidad. Las geometrías las crearemos usando el módulo geom
de la extensión de scripting, en concreto la función createGeometry, que
necesita el tipo de geometría a crear y las dimensiones de la geometría. El
código que podíamos usar es:

geometry_multiline = geom.createGeometry(MULTILINE, D2)

values = dict()
values["ID"] = id #campo calculado equivale a numero de registro
values["LINEA"] = linea_id #primer dato de nuestro rumbo
values["GRADOS"] = grados #grados decimales de nuestro rumbo
values["MINUTOS"] = minutos #minutos decimales de nuestro rumbo
values["SEGUNDOS"] = segundos #segundos decimales de nuestro rumbo
values["DISTANCIA"] = radio #Distancia de nustro rumbo
values["RADIAN"] = angulo #radian calculado
values["X"] = x #Coordenada x calculada
values["Y"] = y #Coordenada y calculada

shape_line.append(values)

Llegados a este punto hemos visto todo lo necesario para obtener lo que
queremos a partir de nuestro rumbo, solo nos falta montar el puzzle y el
resultado sería este.

El resultado final es este

Polígono generado a partir del rumbo

El código final lo podéis descargar desde aquí .

Espero que lo disfrutéis tanto como yo haciéndolo.

¡¡Hasta la próxima!!