12. Módulo de geometrías: geom¶

Funciones mas relevantes que nos encontramos en este modulo:

createGeometry(type, subtype=D2):

Creates an empty geometry

Parámetros:

type (int) – Geometry type
subtype (int) – Geometry dimension

Devuelve:

Geometry

createPoint(subtype=D2, coords=None):

Creates a point geometry

Parámetros:: subtype (int) – Geometry dimension
Devuelve:: Point

createMultiPoint(subtype=D2, points=None):

Creates an multipoint geometry

Parámetros:

subtype (int) – Geometry dimension
points (List of Points) – Points objects

Devuelve:

MultiPoint

createLine(subtype=D2, vertexes=None):

Creates an empty geometry

Parámetros:

subtype (int) – Geometry dimension
vertexes (List of Points)

Devuelve:

Geometry

createMultiLine(subtype=D2, lines=None):

Creates an empty geometry

Parámetros:

subtype (int) – Geometry dimension
lines (List of Lines)

Devuelve:

Geometry

createPolygon(subtype=D2, vertexes=None):

Creates an empty geometry

Parámetros:

subtype (int) – Geometry dimension
vertexes (List of Points)

Devuelve:

Geometry

createMultiPolygon(subtype=D2, polygons=None):

Creates an empty geometry

Parámetros:

subtype (int) – Geometry dimension
polygons (List of Poligons)

Devuelve:

Geometry

createEnvelope(pointMin=None, pointMax=None, dimension=2):

Creates a Envelope 2D

Parámetros:

pointMin (Point) – Geometry
pointMax (Point) – Geometry
dimension (int) – Dimension of the envelope

Devuelve:

Geometry

createGeometryFromWKT(WKT):

Creates an empty geometry

Parámetros:

type (int) – Geometry type
subtype (int) – Geometry dimension

Devuelve:

Geometry

Una geometría es un objeto que contiene información geométrica. Estas geometrías tienen un tipo principal: Point, Line, Polygon.. y un subtipo o dimensión: D2, D3, D2M..

Para el Módulo de Scripting hemos creado la librería gvsig.geom que nos ayudará para crear rápidamente las geometrías que necesitemos. Para algunas operaciones más complicadas tendremos que usar la API de gvSIG desktop.

Para establecer estos tipos y subtipos lo haremos utilizando las constantes que se incluyen en la librería.

Constantes que aparecen en la librería gvsig.geom para la creación de geometrías:

#GeometryTypes
AGGREGATE = Geometry.TYPES.AGGREGATE
ARC = Geometry.TYPES.ARC
CIRCLE = Geometry.TYPES.CIRCLE
CURVE = Geometry.TYPES.CURVE
ELLIPSE = Geometry.TYPES.ELLIPSE
ELLIPTICARC = Geometry.TYPES.ELLIPTICARC
GEOMETRY = Geometry.TYPES.GEOMETRY
MULTICURVE = Geometry.TYPES.MULTICURVE
MULTIPOINT = Geometry.TYPES.MULTIPOINT
MULTISOLID = Geometry.TYPES.MULTISOLID
MULTISURFACE = Geometry.TYPES.MULTISURFACE
NULL = Geometry.TYPES.NULL
POINT = Geometry.TYPES.POINT
SOLID =  Geometry.TYPES.SOLID
SPLINE = Geometry.TYPES.SPLINE
SURFACE = Geometry.TYPES.SURFACE

# Common named geometry types
POLYGON = Geometry.TYPES.SURFACE
LINE = Geometry.TYPES.CURVE
MULTILINE = Geometry.TYPES.MULTICURVE
MULTIPOLYGON = Geometry.TYPES.MULTISURFACE

# geometrySubTypes
D2 = Geometry.SUBTYPES.GEOM2D
D2M = Geometry.SUBTYPES.GEOM2DM
D3 = Geometry.SUBTYPES.GEOM3D
D3M = Geometry.SUBTYPES.GEOM3DM
UNKNOWN = Geometry.SUBTYPES.UNKNOWN

# Dimensions
DIMENSIONS = Geometry.DIMENSIONS

Ejemplo testeando la librería de geom:

import gvsig
reload(gvsig)
from gvsig import *
from gvsig import geom


def main(*args):

    # Create Polygon
    print "\nCreate Polygon"
    x = geom.createPolygon()
    pol_1 = geom.createPolygon(vertexes=[geom.createPoint(geom.D2, 4,5),geom.createPoint(geom.D2,3,3),geom.createPoint(geom.D2,3,2),geom.createPoint(geom.D2,4,5)])
    print "pol_1: ", pol_1.convertToWKT()
    pol_2 = geom.createPolygon(vertexes=[geom.createPoint(geom.D2,4,5),geom.createPoint(geom.D2,3,3),geom.createPoint(geom.D2,3,2),geom.createPoint(geom.D2,4,5)])
    print "pol_2: ", pol_2.convertToWKT()
    pol_3 = geom.createPolygon(vertexes=[geom.createPoint(geom.D2,4,5),geom.createPoint(geom.D2,3,3),geom.createPoint(geom.D2,3,2),geom.createPoint(geom.D2,4,5)])
    print "pol_3: ", pol_3.convertToWKT()
    pol_4 = geom.createPolygon(geom.D2,[(0,0),(10,10),[3,3],[3,6],[0,0]])
    print "pol_4: ", pol_4.convertToWKT()


    # Create MultiPolygon
    print "\nCreate MultiPolygon"
    multipolygon1 = geom.createMultiPolygon()
    multipolygon1.addSurface(pol_1)
    multipolygon1.addSurface(pol_2)
    multipolygon1.addSurface(pol_3)
    print "multipolygon1: ", multipolygon1.convertToWKT()

    multipolygon2 = geom.createMultiPolygon(polygons=[pol_1, pol_2, pol_3])
    print "multipolygon2: ", multipolygon2.convertToWKT()

    p2 = geom.createPoint(geom.D2,1,2)
    print "p2:", p2
    line2 = geom.createLine(geom.D2, [[10,19],p2,[5,2]])
    print line2
    print line2.convertToWKT()

    y = geom.createPoint(geom.D3M,10,1,5,8)
    z = geom.createPoint(geom.D3,10, 1, 5)
    print "point y: ", y,type(y)
    print "point z", z, type(z)

    print "preparing 3d"
    x = geom.createLine(geom.D3M,[(10,10,100,8),(1,95,2,8)])
    print x.convertToWKT()

    # Create point
    print "\nCreate Point"
    point1 = geom.createPoint(geom.D2,10, 10)
    point2 = geom.createGeometry(geom.POINT)
    point2.setX(15)
    point2.setY(15)
    print "Point1: ", point1
    print "Point2: ", point2

    point1 = geom.createPoint(geom.D2,10, 10)
    # Create line
    print "\nCreate Line"
    line1 = geom.createGeometry(geom.LINE)
    line1.addVertex(geom.createPoint(geom.D2,0,0))
    line1.addVertex(geom.createPoint(geom.D2,10,10))
    print "Line1: ", line1.convertToWKT()

    p2 = geom.createPoint(geom.D2,1,2)
    print " === LINE == "
    line2 = geom.createLine(geom.D2, [[10,19], point1 ,p2,[5,2]])
    print "Line2 object: ", line2
    print "Line2: ", line2.convertToWKT()
    print "1", line2.getVertex(0)
    print "2", line2.getVertex(1)
    print "3", line2.getVertex(2)
    print "4", line2.getVertex(3)

    # Create polygon
    print "\nCreate Polygon"
    g = geom.createGeometry(geom.POLYGON)
    g.addVertex(geom.createPoint(geom.D2,0,0))
    g.addVertex(geom.createPoint(geom.D2,10,10))
    g.addVertex(geom.createPoint(geom.D2,10,0))
    g.addVertex(geom.createPoint(geom.D2,0, 0))


    print "JTS of the Polygon 1: ", g.convertToWKT()
    g.setVertex(2, geom.createPoint(geom.D2, 15, 15))
    print "JTS of the Polygon 1 modified: ", g.convertToWKT()

    poli_1 = geom.createPolygon(geom.D2, [[0,0],[1,1],[2,3],[3,6],[0,0]])
    print "Poli_1", poli_1.convertToWKT()
    poli_2 = geom.createPolygon(geom.D3, [[0,1,2],[1,1,5],geom.createPoint(geom.D3,2,1,5),[0,1,2]])
    print "Poli_2", poli_2

    # Create gvSIG geometry from a WKT or WKB
    print "\nCreate gvSIG geometry from WKT or WKB"
    wkt = "POLYGON ((0 0, 150 150, 100 0, 0 0))"
    x = geom.createGeometryFromWKT(wkt)
    print "Polygon from WKT: ", x
    print "Type polygon: ", type(x)

    # Create 3D geometry
    print "\nCreate 3D Geometry"
    p3d = geom.createGeometry(geom.POINT, geom.D3)
    p3d.setX(10)
    p3d.setY(10)
    p3d.setZ(100)
    print "Point 3D: ", p3d, type(p3d)
    p1_3d = geom.createPoint(geom.D3,1,3,3)
    print "Point 3D P1: ", p1_3d

    # Create Multipoint
    print "\nCreate Multipoint: "
    multipoint1 = geom.createMultiPoint(points=[geom.createPoint(geom.D2, 10,10), geom.createPoint(geom.D2,5,2), geom.createPoint(geom.D2,8,3)])
    print "multipoint1: ", multipoint1.convertToWKT()
    multipoint1.addPrimitive(geom.createPoint(geom.D2,3, 2))
    print "multipoint1 modified: ", multipoint1.convertToWKT()

    multipoint2 = geom.createMultiPoint()
    print "multipoint2: ", multipoint2.convertToWKT()

    multipoint3 = geom.createMultiPoint(geom.D3,[[19,10,8],[3,5,7],[35,5,5]])
    print "multipoint3: ", multipoint3.convertToWKT()

    # Create Polygon
    print "\nCreate Polygon"
    x = geom.createPolygon()
    pol_1 = geom.createPolygon(vertexes=[geom.createPoint(geom.D2,4,5),geom.createPoint(geom.D2,3,3),geom.createPoint(geom.D2,3,2),geom.createPoint(geom.D2,4,5)])
    print "pol_1: ", pol_1.convertToWKT()
    pol_2 = geom.createPolygon(vertexes=[geom.createPoint(geom.D2,4,5),geom.createPoint(geom.D2,3,3),geom.createPoint(geom.D2,3,2),geom.createPoint(geom.D2,4,5)])
    print "pol_2: ", pol_2.convertToWKT()
    pol_3 = geom.createPolygon(vertexes=[geom.createPoint(geom.D2,4,5),geom.createPoint(geom.D2,3,3),geom.createPoint(geom.D2,3,2),geom.createPoint(geom.D2,4,5)])
    print "pol_3: ", pol_3.convertToWKT()
    pol_4 = geom.createPolygon(geom.D2,[(0,0),(10,10),[3,3],[3,6],[0,0]])
    print "pol_4: ", pol_4.convertToWKT()


    # Create MultiPolygon
    print "\nCreate MultiPolygon"
    multipolygon1 = geom.createMultiPolygon()
    multipolygon1.addSurface(pol_1)
    multipolygon1.addSurface(pol_2)
    multipolygon1.addSurface(pol_3)
    print "multipolygon1: ", multipolygon1.convertToWKT()

    multipolygon2 = geom.createMultiPolygon(polygons=[pol_1, pol_2, pol_3])
    print "multipolygon2: ", multipolygon2.convertToWKT()

    multipolygon3 = geom.createMultiPolygon(geom.D2,[[[0,0],[1,1],[2,2],[0,0]],[[2,5],[3,5],[1,2],[2,5]],pol_4])
    print "multipolygon3: ", multipolygon3.convertToWKT()

    # CreateLine
    print "\nCreate Line"
    line1 = geom.createLine()


    line1.addVertex(geom.createPoint2D(1,1))
    line1.addVertex(geom.createPoint2D(3,3))
    print "line1: ", line1.convertToWKT()
    line2 = geom.createLine(vertexes=[geom.createPoint2D(0,0), geom.createPoint2D(10,10)])
    print "line2: ", line2.convertToWKT()
    line3 = geom.createLine(geom.D2,[[0,1],[1,5],[5,3]])
    print "line3: ", line3.convertToWKT()

    # Create MultiLine
    print "\nCreate MultiLine"
    multiline1 = geom.createMultiLine()
    multiline1.addCurve(line1)
    multiline1.addCurve(line2)
    print "multiline1: ", multiline1.convertToWKT()

    multiline2 = geom.createMultiLine(lines=[line1, line2], subtype=geom.D2)
    print "multiline2: ", multiline2.convertToWKT()


    # Create Envelope
    envelope = geom.createEnvelope(pointMin=geom.createPoint2D(10,0),pointMax=geom.createPoint2D(10,20))
    print "envelope: ", envelope
    env1 = geom.createEnvelope(point1,[38,29])
    print "env1: ", env1,type(env1)

    # Create from WKT
    print "\nCreate geometr from WKT"
    wkt = geom.createGeometryFromWKT("MULTIPOLYGON (((4 5, 3 3, 3 2, 4 5)), ((4 5, 3 3, 3 2, 4 5)), ((4 5, 3 3, 3 2, 4 5)))")
    print "wkt: ", wkt.convertToWKT()

    # Create 2D
    pg = geom.createPoint2D(19,5)
    print pg
    lg = geom.createLine2D()
    print lg
    mg = geom.createPolygon2D()
    print mg
    mg2 = geom.createPolygon2D([[1,5],[4,5],[1,3],[25,2],pg,[1,5]])
    print mg2.convertToWKT()
    x = geom.createPoint2D(1)
    print x

12.1. Punto¶

Creando puntos:

from gvsig.geom import *

def main(*args):

  print "\nCreate Point"
  point1 = createPoint(D3, 10, 10, 5)
  point2 = createPoint(D2, 15, 7)

  point3 = createGeometry(POINT)
  point3.setX(15)
  point3.setY(15)

  point4 = createPoint(D3M, 5, 18, 3, 2)

  print "Point1: ", point1
  print "Point2: ", point2
  print "Point3: ", point3
  print "Point4: ", point4

Consola:

Create Point
Point1:  POINT Z (10.0 10.0 5.0)
Point2:  POINT (15.0 7.0)
Point3:  POINT (15.0 15.0)
Point4:  POINT ZM (5.0 18.0 3.0 2.0)

Creando punto con 3 dimensiones:

from gvsig.geom import *

def main(*args):

  #Create 3D Geometry
  p3d = createGeometry(POINT, D3)
  p3d.setX(10)
  p3d.setY(10)
  p3d.setZ(100)

  print p3d.convertToWKT()

Consola:

POINT Z (10 10 100)

12.2. Línea¶

Creando líneas:

from gvsig.geom import *

def main(*args):

  line1 = createGeometry(LINE)
  line1.addVertex(createPoint(D2,0,0))
  line1.addVertex(createPoint(D2,10,10))
  print "Line1: ", line1.convertToWKT()

Consola:

Line1:  LINESTRING (0 0, 10 10)

Accediendo a los vértices de la línea:

from gvsig.geom import *

def main(*args):

  line1 = createGeometry(LINE)
  line1.addVertex(createPoint(D2,0,0))
  line1.addVertex(createPoint(D2,10,10))
  print "Line1: ", line1.convertToWKT()

  p1 = createPoint(D2, 3, 01)
  p2 = createPoint(D2, 1, 2)

  # List of vertex or coordinates
  vertx = [[10,19], p1 ,p2,[5,2]]
  line2 = createLine(D2, vertx)

  print "\nLine2 object: ", line2
  print "Line2: ", line2.convertToWKT()
  print "1", line2.getVertex(0)
  print "2", line2.getVertex(1)
  print "3", line2.getVertex(2)
  print "4", line2.getVertex(3)

Consola:

Line1:  LINESTRING (0 0, 10 10)

Line2 object:  Line:2D
Line2:  LINESTRING (10 19, 3 1, 1 2, 5 2)
1 POINT (10.0 19.0)
2 POINT (3.0 1.0)
3 POINT (1.0 2.0)
4 POINT (5.0 2.0)

12.3. Polígono¶

Creando polígonos:

from gvsig.geom import *

def main(*args):

  g = createGeometry(POLYGON, D2)
  g.addVertex(createPoint(D2,0,0))
  g.addVertex(createPoint(D2,10,10))
  g.addVertex(createPoint(D2,10,0))
  g.addVertex(createPoint(D2,0, 0))

  print "WKT Polygon: ", g.convertToWKT()

  poli_1 = createPolygon(D2, [[0,0],[1,1],[2,3],[3,6],[0,0]])
  print "Poli_1", poli_1.convertToWKT()

  poli_2 = createPolygon(D3, [[0,1,2],[1,1,5],createPoint(D3,2,1,5),[0,1,2]])
  print "Poli_2", poli_2.convertToWKT()

Consola:

WKT Polygon:  POLYGON ((0 0, 10 10, 10 0, 0 0))
Poli_1 POLYGON ((0 0, 1 1, 2 3, 3 6, 0 0))
Poli_2 POLYGON Z ((0 1 2, 1 1 5, 2 1 5, 0 1 2))

12.4. Multipunto¶

Creando multipunto:

from gvsig.geom import *

def main(*args):

  # Create Multipoint
  multipoint1 = createMultiPoint(D2, [createPoint(D2, 10,10), createPoint(D2,5,2), createPoint(D2,8,3)])

  print "multipoint1: ", multipoint1.convertToWKT()

  multipoint2 = createMultiPoint()
  print "multipoint2: ", multipoint2.convertToWKT()

  multipoint3 = createMultiPoint(D3,[[19,10,8],[3,5,7],[35,5,5]])
  print "multipoint3: ", multipoint3.convertToWKT()

Consola:

multipoint1:  MULTIPOINT (10 10, 5 2, 8 3)
multipoint2:  MULTIPOINT EMPTY
multipoint3:  MULTIPOINT Z (19 10 8, 3 5 7, 35 5 5)

Añadiendo punto a una geometría de tipo multipunto:

from gvsig.geom import *

def main(*args):

  # Create Multipoint
  multipoint1 = createMultiPoint(D2, [createPoint(D2, 10,10), createPoint(D2,5,2), createPoint(D2,8,3)])

  print "multipoint1: ", multipoint1.convertToWKT()

  multipoint1.addPrimitive(createPoint(D2,3, 2))
  print "multipoint1 modified: ", multipoint1.convertToWKT()

Consola:

multipoint1:  MULTIPOINT (10 10, 5 2, 8 3)
multipoint1 modified:  MULTIPOINT (10 10, 5 2, 8 3, 3 2)

12.5. Multilínea¶

Creando multilíneas:

from gvsig.geom import *

def main(*args):

# CreateLine
  print "\nCreate Line"
  line1 = createLine()

  line1.addVertex(createPoint2D(1,1))
  line1.addVertex(createPoint2D(3,3))
  line2 = createLine(vertexes=[createPoint2D(0,0), createPoint2D(10,10)])
  line3 = createLine(D2,[[0,1],[1,5],[5,3]])

  # Create MultiLine
  multiline1 = createMultiLine()
  multiline1.addCurve(line1)
  multiline1.addCurve(line2)
  print "multiline1: ", multiline1.convertToWKT()

  multiline2 = createMultiLine(D2, lines=[line1, line2])
  print "multiline2: ", multiline2.convertToWKT()

  multiline3 = createGeometry(MULTICURVE, D2)
  multiline3.addCurve(line1)
  multiline3.addCurve(line3)
  print "multiline3: ", multiline3.convertToWKT()

Consola:

Create Line
multiline1:  MULTILINESTRING ((1 1, 3 3), (0 0, 10 10))
multiline2:  MULTILINESTRING ((1 1, 3 3), (0 0, 10 10))
multiline3:  MULTILINESTRING ((1 1, 3 3), (0 1, 1 5, 5 3))

12.6. Multipolígono¶

Creando multipolígonos:

from gvsig.geom import *

def main(*args):

  # Create Polygon
  x = createPolygon()
  pol_1 = createPolygon(vertexes=[createPoint(D2,4,5),createPoint(D2,3,3),createPoint(D2,3,2),createPoint(D2,4,5)])
  pol_2 = createPolygon(vertexes=[createPoint(D2,4,5),createPoint(D2,3,3),createPoint(D2,3,2),createPoint(D2,4,5)])
  pol_3 = createPolygon(vertexes=[createPoint(D2,4,5),createPoint(D2,3,3),createPoint(D2,3,2),createPoint(D2,4,5)])
  pol_4 = createPolygon(D2,[(0,0),(10,10),[3,3],[3,6],[0,0]])


  # Create MultiPolygon
  multipolygon1 = createMultiPolygon()
  multipolygon1.addSurface(pol_1)
  multipolygon1.addSurface(pol_2)
  multipolygon1.addSurface(pol_3)
  print "multipolygon1: ", multipolygon1.convertToWKT()

  multipolygon2 = createMultiPolygon(polygons=[pol_1, pol_2, pol_3]) #2D as default
  print "multipolygon2: ", multipolygon2.convertToWKT()

  multipolygon3 = createMultiPolygon(D2, [[[0,0],[1,1],[2,2],[0,0]],
                      [[2,5],[3,5],[1,2],[2,5]],
                      pol_4])
  print "multipolygon3: ", multipolygon3.convertToWKT()

  multipolygon4 = createMultiPolygon()
  print "multipolygon4: ", multipolygon4.convertToWKT()

Consola:

multipolygon1:  MULTIPOLYGON (((4 5, 3 3, 3 2, 4 5)), ((4 5, 3 3, 3 2, 4 5)), ((4 5, 3 3, 3 2, 4 5)))
multipolygon2:  MULTIPOLYGON (((4 5, 3 3, 3 2, 4 5)), ((4 5, 3 3, 3 2, 4 5)), ((4 5, 3 3, 3 2, 4 5)))
multipolygon3:  MULTIPOLYGON (((0 0, 1 1, 2 2, 0 0)), ((2 5, 3 5, 1 2, 2 5)), ((0 0, 10 10, 3 3, 3 6, 0 0)))
multipolygon4:  MULTIPOLYGON EMPTY

12.7. Envelope¶

Crear Envelope del tipo Envelope2D

from gvsig.geom import *

def main(*args):

  # Create Envelope
  envelope = createEnvelope(pointMin=createPoint2D(10,0),pointMax=createPoint2D(10,20))
  print "envelope: ", envelope
  point1 = createPoint(D2, 0, 0)
  env1 = createEnvelope(point1,[38,29])
  print "env1: ", env1
  print "env1 type: ", type(env1)

Consola:

envelope:  POLYGON ((10.0 0.0, 10.0 20.0, 10.0 20.0, 10.0 0.0, 10.0 0.0))
env1:  POLYGON ((0.0 0.0, 0.0 29.0, 38.0 29.0, 38.0 0.0, 0.0 0.0))
env1 type:  <type 'org.gvsig.fmap.geom.jts.primitive.Envelope2D'>

12.8. WKT¶

Crear geometría desde un WKT:

from gvsig.geom import *

def main(*args):

  # Create from WKT
  wkt = createGeometryFromWKT("MULTIPOLYGON (((4 5, 3 3, 3 2, 4 5)), ((4 5, 3 3, 3 2, 4 5)), ((4 5, 3 3, 3 2, 4 5)))")
  print "wkt: ", wkt.convertToWKT()

Consola:

wkt:  MULTIPOLYGON (((4 5, 3 3, 3 2, 4 5)), ((4 5, 3 3, 3 2, 4 5)), ((4 5, 3 3, 3 2, 4 5)))

12.9. Operaciones espaciales¶

Puedes consultar todas las operaciones espaciales en el interfaz de geometrías Geometry

Distancia entre puntos:

from gvsig.geom import *

def main(*args):
  # Creamos punto2
  point1 = createPoint(D2, 0, 0)
  point2 = createPoint(D2, 10, 10)
  print "Distance 2D: ", point1.distance(point2)

  point3 = createPoint(D3, 0, 0, 100)
  point4 = createPoint(D3, 10, 10, 0)
  print "Distance 2D: ", point3.distance(point4)

Consola:

Distance 2D:  14.1421356237
Distance 2D:  14.1421356237

Moviendo un punto:

from gvsig.geom import *

def main(*args):
  # Creamos punto2
  point1 = createPoint(D2, 10, 10)

  print "Point 1: ", point1.convertToWKT()

  #Move point
  point1.move(5, -3)

  print "Moved point by 5, -3: ", point1.convertToWKT()

Consola:

Point 1:  POINT (10 10)
Moved point by 5, -3:  POINT (15 7)

Operaciones entre polígonos y líneas:

from gvsig.geom import *

def main(*args):
  # Creamos punto2
  point1 = createPoint(D2, 0, 0)
  buffer1 = point1.buffer(10)

  line1 = createLine(D2, [[-5, -5],[10, 10]])

  print "Intersects?: ", buffer1.intersects(line1)
  print "Intersection: ", buffer1.intersection(line1).convertToWKT()

Consola:

Intersects?:  True
Intersection:  LINESTRING (-5 -5, 7.071067811865473 7.071067811865473)

Operaciones espaciales entre polígonos:

from gvsig.geom import *

def main(*args):
  # Creamos punto2
  point1 = createPoint(D2, 0, 0)

  # Aplicamos un area de influencia buffer(m)
  buffer1 = point1.buffer(10)
  print "Buffer1 Area: ", buffer1.area()

  # Creamos punto 2
  point2 = createPoint(D2, 8, 0)
  buffer2 = point2.buffer(5)
  print "Buffer2 Area: ", buffer2.area()

  # Union
  buffer12_union = buffer1.union(buffer2)
  print "Buffer12 Union Area: ", buffer12_union.area()

  # Differencia
  buffer12_diff = buffer1.difference(buffer2)
  print "Buffer12 Difference Area: ", buffer12_diff.area()

Consola:

Buffer1 Area:  312.144515226
Buffer2 Area:  78.0361288065
Buffer12 Union Area:  335.886462528
Buffer12 Difference Area:  257.850333721

Creando un polígono, aplicarle un área de influencia, añadirle un anillo interno y agregarlo en una capa nueva:

from gvsig import *
from gvsig.geom import *

def main(*args):

  pol_1 = createPolygon(D2M, [(0,0),(300,0),(300,300),(0,300),(0,0)])


  # Add interior ring
  pol_1x = createPolygon(D2M, pol_1).buffer(200)
  pol_1x.addInteriorRing(pol_1)

  schema = createSchema()
  schema.append("ID", "STRING", 5)
  schema.append("GEOMETRY", "GEOMETRY")
  schema.get('GEOMETRY').setGeometryType(POLYGON, D2M)

  shape = createShape(schema ,CRS='EPSG:25830')

  shape.edit()
  shape.append(ID=1, GEOMETRY=pol_1x)
  shape.commit()

  currentView().addLayer(shape)

Ejemplo para extraer todos los vértices de una capa de polígonos en forma de multipuntos:

# encoding: utf-8

from gvsig import *
from gvsig import geom

def main(*args):
    """ Extraer vertices de poligonos como multipuntos """

    layer = currentLayer()
    features = layer.features()

    sch = createFeatureType()
    sch.append("GEOMETRY", "GEOMETRY")
    sch.get("GEOMETRY").setGeometryType(geom.MULTIPOINT, geom.D2)
    shp = createShape(sch)

    for feature in features:
        gf = feature.geometry()
        shp.append(GEOMETRY=gf.toPoints())

    shp.commit()
    currentView().addLayer(shp)

Resultado:

Con una pequeña modificación, podemos hacer que añada un punto por cada vértice, recorriendo las geometrías de multipunto y extrayendo uno a uno sus puntos:

# encoding: utf-8

from gvsig import *
from gvsig import geom

def main(*args):
    """ Extraer vertices de poligonos como multipuntos """

    layer = currentLayer()
    features = layer.features()

    sch = createFeatureType()
    sch.append("GEOMETRY", "GEOMETRY")
    sch.get("GEOMETRY").setGeometryType(geom.POINT, geom.D2)
    shp = createShape(sch)

    for feature in features:
        gf = feature.geometry()

        for p in gf:
            shp.append(GEOMETRY=p)

    shp.commit()
    currentView().addLayer(shp)

Otra forma de recorrer los vertices podría ser haciendo uso del método getVertex y getNumVertices:

# encoding: utf-8

from gvsig import *
from gvsig import geom

def main(*args):
    """ Extraer vertices de poligonos como multipuntos """

    layer = currentLayer()
    features = layer.features()

    for feature in features:
        gf = feature.geometry()
        numvertices = gf.getNumVertices()
        print "Vertices: ", numvertices
        for i in range(0, numvertices):
            print gf.getVertex(i) #geometria de tipo punto

Dando por consola una salida similar a:

Vertices:  7
POINT (274.17769116287 194.9280163662371)
POINT (275.81171697787 166.3325646037371)
POINT (248.03327812287003 105.05659654123708)
POINT (145.90666468537003 40.5125768487371)
POINT (162.24692283537004 128.74997085873707)
POINT (207.99964565537007 140.1881515637371)
POINT (274.17769116287 194.9280163662371)
...

Aquí un ejemplo de creación de geometrías por su dimensión:

# encoding: utf-8

import gvsig
from gvsig import geom
from org.gvsig.fmap.geom import Geometry
from org.gvsig.fmap.geom import GeometryLocator

def main(*args):

    geomManager = GeometryLocator.getGeometryManager()
    print "\nCreate Point"
    point1 = geomManager.create(Geometry.TYPES.POINT, Geometry.SUBTYPES.GEOM3DM)
    print point1
    dimension = point1.getDimension()
    for i in range(0,point1.getDimension()):
        point1.setCoordinateAt(i,10)
    print point1

Un ejemplo avanzado de recorrer todos los vértices de una geometría, sea cual sea tu tipo, anillos que contiene, o dimensión de cada una, para copiarlo a otra geometría vacía a la vez que se modifica su dimensión a D2:

# encoding: utf-8

import gvsig
from gvsig import geom
from org.gvsig.fmap.mapcontext.layers.vectorial import SpatialEvaluatorsFactory

from org.gvsig.expressionevaluator import ExpressionEvaluatorLocator
from org.gvsig.fmap.geom.aggregate import MultiPrimitive
from org.gvsig.fmap.geom.primitive import Polygon, Point
from org.gvsig.fmap.geom import GeometryLocator
from org.gvsig.fmap.geom.aggregate import MultiPrimitive

from org.gvsig.fmap.geom import Geometry

def transformPointTo2D(iVertex,nv=None):
    if nv==None:
        nv = GeometryLocator.getGeometryManager().create(iVertex.getGeometryType().getType(),Geometry.SUBTYPES.GEOM2D)
    for d in range(0,iVertex.getDimension()):
        try:
            nv.setCoordinateAt(d,iVertex.getCoordinateAt(d))
        except:
            pass
    return nv

def insertVertexFromGeometryInGeometry(iPol,np,transformMethod=None):
    geomManager = GeometryLocator.getGeometryManager()
    if isinstance(iPol, Point):
        if transformMethod==None:
            for d in range(0,iPol.getDimension()):
                np.setCoordinateAt(d,iPol.getCoordinateAt(d))
            return
        else:
            transformMethod(iPol,np)
            return

    for v in range(0, iPol.getNumVertices()):
        iVertex = iPol.getVertex(v)
        if transformMethod==None:
            nv = geomManager.create(iVertex.getGeometryType().getType(),iVertex.getGeometryType().getSubType())
            for d in range(0,iVertex.getDimension()):
                nv.setCoordinateAt(d,iVertex.getCoordinateAt(d))
        else:
            nv = transformMethod(iVertex)
        np.addVertex(nv)

    if isinstance(iPol, Polygon):
        for r in range(0, iPol.getNumInteriorRings()):
            iRing = iPol.getInteriorRing(r)
            nr = geomManager.create(iRing.getGeometryType().getType(),iRing.getGeometryType().getSubType())
            insertVertexFromGeometryInGeometry(iRing, nr,transformMethod)
            np.addInteriorRing(nr)

def process():
    store = gvsig.currentLayer().getFeatureStore()
    ## method of transformation
    method = "transformPointTo2D"
    if method == "transformPointTo2D":
        transformMethod = transformPointTo2D
        subtype = geom.D2
    # process
    geomManager = GeometryLocator.getGeometryManager()
    fset = store.getFeatureSet()
    nsch = gvsig.createFeatureType(store.getDefaultFeatureType())
    subtype = nsch.get("GEOMETRY").getGeometrySubType()
    nsch.get("GEOMETRY").setGeometryType(nsch.get("GEOMETRY").getGeometryType(), subtype)
    outputFilePath = gvsig.getTempFile("result_geometries",".shp")

    ns = gvsig.createShape(nsch,outputFilePath)
    ns.edit()
    store = ns.getFeatureStore()
    for f in fset:
        fg = f.getDefaultGeometry()
        if subtype == None:
            subtype =  fg.getGeometryType().getSubType()
        newGeometry = geomManager.create(fg.getGeometryType().getType(), subtype)
        if isinstance(fg,MultiPrimitive):
            for i in range(0,fg.getPrimitivesNumber()):
                iPol = fg.getPrimitiveAt(i)
                np = geomManager.create(iPol.getGeometryType().getType(), subtype)
                insertVertexFromGeometryInGeometry(iPol, np, transformMethod)
                newGeometry.addPrimitive(np)
        else:
            insertVertexFromGeometryInGeometry(fg, newGeometry,transformMethod)
        nf = store.createNewFeature(f)
        nf.set("GEOMETRY", newGeometry)
        store.insert(nf)
    ns.finishEditing()
    gvsig.currentView().addLayer(ns)

def main(*args):
    process()

12.10. Transformación entre proyecciones¶

En el siguiente ejemplo vemos como transformar una geometría entre dos proyecciones:

# encoding: utf-8

from gvsig import *
from gvsig.geom import *

def main(*args):
    """ Convertir geometria entre proyecciones """

    crs1 = getCRS('EPSG:4326')
    crs2 = getCRS('EPSG:32630')

    #Get IProjection
    #view1 = currentProject().getViews()[0].getProjection()
    #view2 = currentProject().getViews()[1].getProjection()

    #Get ICoordTrans
    ICoordTrans1 = crs1.getCT(crs2)

    point = createPoint2D(-0.375,39.466) #Valencia
    print "point: ", point
    point.reProject(ICoordTrans1)
    print "reprojected: ", point

12.11. isSubtype¶

Si queremos identificar si una capa es de un tipo o de otro, podemos hacerlo utilizando el GeometryManager.

Para ello, se busca el tipo de geometría que tiene una capa y se le pregunta si es instancia de alguna de las que queremos.

En este caso, que sean geometrías de tipo SURFACE:

# encoding: utf-8

import gvsig
from org.gvsig.fmap.geom import Geometry
from org.gvsig.fmap.geom import GeometryLocator

# Con geometrias normales se quedaria con el getGeometryType()
def main(*args):
view = gvsig.currentView()
geomManager = GeometryLocator.getGeometryManager()
for layer in view.deepiterator():
    getTypeVectorLayer = getattr(layer,"getTypeVectorLayer",None)
    if getTypeVectorLayer != None:
    geomType = getTypeVectorLayer()
    print geomType.getName()
    # depende lo fino que quieras hilar puedes sustituir SURFACE por POLYGON
    if (geomManager.isSubtype(Geometry.TYPES.SURFACE,geomType.getType()) or
        geomManager.isSubtype(Geometry.TYPES.MULTISURFACE,geomType.getType())):
        print layer

12. Módulo de geometrías: geom¶

12.1. Punto¶

12.2. Línea¶

12.3. Polígono¶

12.4. Multipunto¶

12.5. Multilínea¶

12.6. Multipolígono¶

12.7. Envelope¶

12.8. WKT¶

12.9. Operaciones espaciales¶

12.10. Transformación entre proyecciones¶

12.11. isSubtype¶

Tabla de contenido

Esta página

Links