gvSIG blog

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:

Una única entidad seleccionada.

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

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.