Partir polilíneas y dividirlas por cada intersección mediante Scripting en gvSIG

Un geoproceso que podemos necesitar es el de partir polilíneas en líneas simples, y además, generar una división en ellas por cada intersección con otra línea. Este geoproceso se ejecuta sobre toda la capa, comprobando las posibles intersecciones entre todas las polilíneas. Lo podemos necesitar para calcular intersecciones, longitudes, selecciones, etc.

Ejecutaríamos el script teniendo un fichero lleno de polilíneas que se superponen similar al de la imagen:

2015-06-29 14_28_40-gvSIG 2.2.0.2308 RC1 _ Sin título

Una única entidad seleccionada.

Obtendríamos un resultado como el siguiente, una capa de líneas partidas:

2015-06-29 14_32_30-gvSIG 2.2.0.2308 RC1 _ Sin título

Varias entidades seleccionadas divisiones de la anterior. Todas las polilíneas están partidas de igual modo.

Este geoproceso es similar al correspondiente utilizado por GRASS de ‘v.clean-break’.

Si estás interesado en aprender a realizar estas operaciones espaciales y mucho más, apúntate ya en el MOOC de Scripting en gvSIG 2. Curso gratuito y en abierto, con mucha actividad en los foros, siendo opcional la obtención del certificado.

El código del script es el siguiente, el único valor que deberás de cambiar es el correspondiente a la ruta del fichero de salida.

import gvsig
import os
from geom import *

def main():
    clean_break(gvsig.currentLayer(), "C:/temp/")

def clean_break(layer, path):
    output = crearCapa_cleanBreak(layer, path)

    features = layer.features()
    for f in features:
        geom = f.geometry()
        values = f.getValues()
        numVertex = geom.getNumVertices()

        for i in range(0,int(numVertex) - 1):
            ver1 = geom.getVertex(i)
            ver2 = geom.getVertex(i+1)
            listVertex = []
            listVertex.append((ver1, 0))
            listVertex.append((ver2, ver2.distance(ver1)))
            lineGoingDivided = createGeometry(2)
            lineGoingDivided.addVertex(ver1)
            lineGoingDivided.addVertex(ver2)
            for n in layer.features(): #getSelection():#.features():
                    geomN = n.geometry()
                    numVertexN = geomN.getNumVertices()
                    for ni in range(0,int(numVertexN) - 1):
                        ver1n = geomN.getVertex(ni)
                        ver2n = geomN.getVertex(ni+1)
                        lineGoingCross = createGeometry(2)
                        lineGoingCross.addVertex(ver1n)
                        lineGoingCross.addVertex(ver2n)
                        if lineGoingCross == lineGoingDivided:
                            pass
                        if lineGoingCross.crosses(lineGoingDivided):
                            point = lineGoingCross.intersection(lineGoingDivided)
                            if str(point)== "MultiPoint2D":
                                print "Error: Interseccion multipunto"
                                #nGeometries = point.getGeometries()
                                #for nG in nGeometries:
                                #    if lineGoingDivided.crosses(nG) and n.crosses(nG):
                                #        listVertex.append((nG, nG.distance(ver1)))
                            elif point.type == 1: #de tipo punto
                                #print "Inserción normal"
                                listVertex.append((point, point.distance(ver1)))
                            else:
                                print "Error Pass, geometria tipo: ", point.type
                                pass

            sortedVertex = list(sorted(list(set(listVertex)), key=lambda vertex: vertex[1]))

            # "Creacion de menores"
            if len(sortedVertex) >= 2:
                for l in range(0, len(sortedVertex)-1):
                    tramo = createGeometry(2)
                    tramo.addVertex(sortedVertex[l][0])
                    tramo.addVertex(sortedVertex[l+1][0])
                    values["GEOMETRY"] = tramo
                    output.append(values)

    output.commit()
    gvsig.currentView().addLayer(output)
    return output

def crearCapa_cleanBreak(layer, path="C:/temp/"):
   nombre = "clean_break"
   numero = 0
   #TIP
   while True:
       archivo = nombre + str(numero)
       ruta = path + archivo + ".shp"
       if os.path.exists(ruta):
           numero +=1
       else: break
   #/TIP

   CRS = gvsig.currentProject().getProjectionCode()
   schema = gvsig.createSchema(layer.getSchema())
   #schema.append("GEOMETRY","GEOMETRY")
   output = gvsig.createShape ( schema, ruta, CRS=CRS,
       geometryType= LINE)
   return output

Puedes ver cómo crear el script en gvSIG en el siguiente vídeo sobre cómo cargar un script desde el repositorio.

This entry was posted in development, gvSIG Desktop, scripting, spanish. Bookmark the permalink.

2 Responses to Partir polilíneas y dividirlas por cada intersección mediante Scripting en gvSIG

  1. Pingback: Partir polilíneas y dividirlas por cada intersección mediante Scripting en gvSIG | Geo-How-To News

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s