Clipping polylines and creating intersections using Scripting on gvSIG

There’s a useful geoprocess to clip simple lines and create a separation at the intersections with other lines. This geoprocess is run on the whole layer, checking the intersections between all the polylines. We can need it to calculate intersection, lengths, selections…

We would run the script when we have a file with polylines similar to this image:

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

Una única entidad seleccionada.

We would get a results like this one, a layer with clipped lines:

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.

This geoprocess is similar to the ‘v.clean-break’ tool of GRASS.

The source code of the script is this one, you only will have to change the path to the output file.

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 "Normal insertion"
                                listVertex.append((point, point.distance(ver1)))
                            else:
                                print "Error Pass, geometry type: ", 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

You can see how to create the script in gvSIG at the next video about how to load a script from the repository.

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

One Response to Clipping polylines and creating intersections using Scripting on gvSIG

  1. Pingback: Clipping polylines and creating intersections using Scripting on #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 )

Twitter picture

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

Facebook photo

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

Google+ photo

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

Connecting to %s