example.py

在Linux下做的QuadTree的程序

import os
import sys
import distutils.util
import math

# Put local build directory on head of python path
platformdir = '-'.join((distutils.util.get_platform(),
                        '.'.join(map(str, sys.version_info[0:2]))))
sys.path.insert(0, os.path.join('build', 'lib.' + platformdir))

# import geos from the local build directory
import geos

pm = geos.PrecisionModel(2.0,0,0)
global_factory = geos.GeometryFactory(pm,-1)

def wkt_print_geoms(geoms):
    wkt = geos.WKTWriter()
    size = len(geoms)
    for i in range(size):
        tmp = wkt.write(geoms[i])
        print "[%d]" % i, tmp

def create_point(x,y):
    c = geos.Coordinate(x,y)
    p = global_factory.createPoint(c)

    return p

def create_ushaped_linestring(xoffset, yoffset, side):
    cl = geos.DefaultCoordinateSequence()

    cl.add(geos.Coordinate(xoffset, yoffset))
    cl.add(geos.Coordinate(xoffset, yoffset+side))
    cl.add(geos.Coordinate(xoffset+side, yoffset+side))
    cl.add(geos.Coordinate(xoffset+side, yoffset))

    ls = global_factory.createLineString(cl)
    return ls

def create_square_linearring(xoffset,yoffset,side):
    cl = geos.DefaultCoordinateSequence()

    cl.add(geos.Coordinate(xoffset,yoffset))
    cl.add(geos.Coordinate(xoffset,yoffset+side))
    cl.add(geos.Coordinate(xoffset+side,yoffset+side))
    cl.add(geos.Coordinate(xoffset+side,yoffset))
    cl.add(geos.Coordinate(xoffset,yoffset))

    lr = global_factory.createLinearRing(cl)
    return lr

def create_square_polygon(xoffset,yoffset,side):
    outer = create_square_linearring(xoffset,yoffset,side)
    inner = create_square_linearring(xoffset+(side/3.),yoffset+(side/3.),(side/3.))

    holes = geos.vector_GeometryP()
    holes.push_back(inner)

    poly = global_factory.createPolygon(outer,holes)
    return poly

def create_simple_collection(geoms):
    collect = geos.vector_GeometryP()
    for i in geoms:
        collect.push_back(i)
    return global_factory.createGeometryCollection(collect)

def create_circle(centerX,centerY,radius):
    shapefactory = geos.GeometricShapeFactory(global_factory)
    shapefactory.setCentre(geos.Coordinate(centerX, centerY))
    shapefactory.setSize(radius)
    return shapefactory.createCircle()

def create_ellipse(centerX,centerY,width,height):
    shapefactory = geos.GeometricShapeFactory(global_factory)
    shapefactory.setCentre(geos.Coordinate(centerX, centerY))
    shapefactory.setHeight(height)
    shapefactory.setWidth(width)
    return shapefactory.createCircle()

def create_rectangle(llX,llY,width,height):
    shapefactory = geos.GeometricShapeFactory(global_factory)
    shapefactory.setBase(geos.Coordinate(llX, llY))
    shapefactory.setHeight(height)
    shapefactory.setWidth(width)
    shapefactory.setNumPoints(4)
    return shapefactory.createRectangle()

def create_arc(llX,llY,width,height,startang,endang):
    shapefactory = geos.GeometricShapeFactory(global_factory)
    shapefactory.setBase(geos.Coordinate(llX, llY))
    shapefactory.setHeight(height)
    shapefactory.setWidth(width)
    #shapefactory.setNumPoints(100) #the default (100 pts)
    return shapefactory.createArc(startang, endang)

def do_all():
    geoms = []

    geoms.append(create_point(150, 350))
    geoms.append(create_ushaped_linestring(60,60,100))
    geoms.append(create_square_linearring(0,0,100))
    geoms.append(create_square_polygon(0,200,300))
    geoms.append(create_square_polygon(0,250,300))
    geoms.append(create_simple_collection(geoms))

    # These ones use a GeometricShapeFactory
    geoms.append(create_circle(0, 0, 10))
    geoms.append(create_ellipse(0, 0, 8, 12))
    geoms.append(create_rectangle(-5, -5, 10, 10)) # a square
    geoms.append(create_rectangle(-5, -5, 10, 20)) # a rectangle

    # The upper-right quarter of a vertical ellipse
    geoms.append(create_arc(0, 0, 10, 20, 0, math.pi/2))

    print "--------HERE ARE THE BASE GEOMS ----------"
    wkt_print_geoms(geoms)


####################
# UNARY OPERATIONS                 #
####################

    ################
    #CENTROID                        #
    ################
    # Find centroid of each base geometry
    newgeoms = []
    for i in range(len(geoms)):
        newgeoms.append(geoms[i].getCentroid())

    print "\n","------- AND HERE ARE THEIR CENTROIDS -----"
    wkt_print_geoms(newgeoms)

    ################
    # BUFFER                           #
    ################
    newgeoms = []
    for i in range(len(geoms)):
        try:
            newgeoms.append(geoms[i].buffer(10))
        except geos.GEOSException():
            exc = geos.GEOSException()
            print "GEOS Exception: geometry ",geoms[i],"->buffer(10): ",exc.toString()

    print "\n","--------HERE COMES THE BUFFERED GEOMS ----------"
    wkt_print_geoms(newgeoms)

    ################
    # CONVEX HULL                 #
    ################
    newgeoms = []
    for i in range(len(geoms)):
        newgeoms.append(geoms[i].convexHull())

    print "\n","--------HERE COMES THE HULLS----------"
    wkt_print_geoms(newgeoms)

####################
# RELATIONAL OPERATORS        #
####################

    print "-------------------------------------------------------------------------------"
    print "RELATIONAL OPERATORS"
    print "-------------------------------------------------------------------------------"

    size = len(geoms)
    ################
    # DISJOINT                          #
    ################

    print
    print "\t".join(["   DISJOINT   "]+["[%d]" % i for i in range(size)])

    for i in range(size):
        print "      [%d]\t" % i,
        for j in range(size):
            try:
                if geoms[i].disjoint(geoms[j]):
                    print " 1\t",
                else:
                    print " 0\t",
            except geos.GEOSException():
                exc = geos.GEOSException()
                print exc.toString()
            except:
                print " X\t",
        print

    ################
    # TOUCHES                         #
    ################

    print
    print "\t".join(["    TOUCHES   "]+["[%d]" % i for i in range(size)])

    for i in range(size):
        print "      [%d]\t" % i,
        for j in range(size):
            try:
                if geoms[i].touches(geoms[j]):
                    print " 1\t",
                else:
                    print " 0\t",
            except geos.GEOSException():
                exc = geos.GEOSException()
                print exc.toString()
            except:
                print " X\t",
        print

    ################
    # INTERSECTS                   #
    ################

    print
    print "\t".join([" INTERSECTS   "]+["[%d]" % i for i in range(size)])

    for i in range(size):
        print "      [%d]\t" % i,
        for j in range(size):
            try:
                if geoms[i].intersects(geoms[j]):
                    print " 1\t",
                else:
                    print " 0\t",
            except geos.GEOSException():
                exc = geos.GEOSException()
                print exc.toString()
            except:
                print " X\t",
        print

    ################
    # CROSSES                        #
    ################

    print
    print "\t".join(["    CROSSES   "]+["[%d]" % i for i in range(size)])

    for i in range(size):
        print "      [%d]\t" % i,
        for j in range(size):
            try:
                if geoms[i].crosses(geoms[j]):
                    print " 1\t",
                else:
                    print " 0\t",
            except geos.GEOSException():
                exc = geos.GEOSException()
                print exc.toString()
            except:
                print " X\t",
        print

    ################
    # WITHIN                        #
    ################

    print
    print "\t".join(["     WITHIN   "]+["[%d]" % i for i in range(size)])

    for i in range(size):
        print "      [%d]\t" % i,
        for j in range(size):
            try:
                if geoms[i].within(geoms[j]):
                    print " 1\t",
                else:
                    print " 0\t",
            except geos.GEOSException():
                exc = geos.GEOSException()
                print exc.toString()
            except:
                print " X\t",
        print

    ################
    # CONTAINS                        #
    ################

    print
    print "\t".join(["   CONTAINS   "]+["[%d]" % i for i in range(size)])

    for i in range(size):
        print "      [%d]\t" % i,
        for j in range(size):
            try:
                if geoms[i].contains(geoms[j]):
                    print " 1\t",
                else:
                    print " 0\t",
            except geos.GEOSException():
                exc = geos.GEOSException()
                print exc.toString()
            except:
                print " X\t",
        print

    ################
    # OVERLAPS                        #
    ################

    print
    print "\t".join(["   OVERLAPS   "]+["[%d]" % i for i in range(size)])

    for i in range(size):
        print "      [%d]\t" % i,
        for j in range(size):
            try:
                if geoms[i].overlaps(geoms[j]):
                    print " 1\t",
                else:
                    print " 0\t",
            except geos.GEOSException():
                exc = geos.GEOSException()
                print exc.toString()
            except:
                print " X\t",
        print

    ################
    # RELATE                            #
    ################

    print
    print "\t".join(["     RELATE   "]+["[%d]" % i for i in range(size)])

    for i in range(size):
        print "      [%d]\t" % i,
        for j in range(size):
            im = geos.IntersectionMatrix('')
            try:
                if geoms[i].relate(geoms[j],"212101212"):
                    print " 1\t",
                else:
                    print " 0\t",
                im=geoms[i].relate(geoms[j])
            except geos.GEOSException():
                exc = geos.GEOSException()
                print exc.toString()
            except:
                print " X\t",
        print

    ################
    # EQUALS                        #
    ################

    print
    print "\t".join(["     EQUALS   "]+["[%d]" % i for i in range(size)])

    for i in range(size):
        print "      [%d]\t" % i,
        for j in range(size):
            try:
                if geoms[i].equals(geoms[j]):
                    print " 1\t",
                else:
                    print " 0\t",
            except geos.GEOSException():
                exc = geos.GEOSException()
                print exc.toString()
            except:
                print " X\t",
        print

    ################
    # EQUALS_EXACT             #
    ################

    print
    print "\t".join(["EQUALS_EXACT  "]+["[%d]" % i for i in range(size)])

    for i in range(size):
        print "      [%d]\t" % i,
        for j in range(size):
            try:
                if geoms[i].equalsExact(geoms[j],0.5):
                    print " 1\t",
                else:
                    print " 0\t",
            except geos.GEOSException():
                exc = geos.GEOSException()
                print exc.toString()
            except:
                print " X\t",
        print

    ################
    # IS_WITHIN_DISTANCE     #
    ################

    print
    print "\t".join(["IS_WITHIN_DIST"]+["[%d]" % i for i in range(size)])

    for i in range(size):
        print "      [%d]\t" % i,
        for j in range(size):
            try:
                if geoms[i].isWithinDistance(geoms[j],2):
                    print " 1\t",
                else:
                    print " 0\t",
            except geos.GEOSException():
                exc = geos.GEOSException()
                print exc.toString()
            except:
                print " X\t",
        print

####################
# COMBINATIONS
####################

    print
    print "-------------------------------------------------------------------------------"
    print "COMBINATIONS"
    print "-------------------------------------------------------------------------------"

    ################
    # UNION
    ################

    newgeoms = []
    for i in range(size-1):
        for j in range(i+1,size):
            try:
                newgeoms.append(geoms[i].Union(geoms[j]))
            except geos.GEOSException():
                exc = geos.GEOSException()
                print exc.toString()
            except:
                pass

    print "\n", "----- AND HERE ARE SOME UNION COMBINATIONS ------"
    wkt_print_geoms(newgeoms)

    ################
    # INTERSECTION
    ################

    newgeoms = []
    for i in range(size-1):
        for j in range(i+1,size):
            try:
                newgeoms.append(geoms[i].intersection(geoms[j]))
            except geos.GEOSException():
                exc = geos.GEOSException()
                print exc.toString()
            except:
                pass

    print "\n", "----- HERE ARE SOME INTERSECTIONS COMBINATIONS ------"
    wkt_print_geoms(newgeoms)

    ################
    # DIFFERENCE
    ################

    newgeoms = []
    for i in range(size-1):
        for j in range(i+1,size):
            try:
                newgeoms.append(geoms[i].difference(geoms[j]))
            except geos.GEOSException():
                exc = geos.GEOSException()
                print exc.toString()
            except:
                pass

    print "\n", "----- HERE ARE SOME DIFFERENCE COMBINATIONS ------"
    wkt_print_geoms(newgeoms)

    ################
    # SYMMETRIC DIFFERENCE
    ################

    newgeoms = []
    for i in range(size-1):
        for j in range(i+1,size):
            try:
                newgeoms.append(geoms[i].symDifference(geoms[j]))
            except geos.GEOSException():
                exc = geos.GEOSException()
                print exc.toString()
            except:
                pass

    print "\n", "----- HERE ARE SYMMETRIC DIFFERENCES ------"
    wkt_print_geoms(newgeoms)

    ################
    # LINEMERGE
    ################
    temp = geos.vector_GeometryP()
    for g in geoms:
        temp.push_back(g)
    lm = geos.LineMerger()
    lm.add(temp)
    mls = lm.getMergedLineStrings()
    newgeoms = []
    for i in range(mls.size()):
        newgeoms.append(mls[i])
    del mls

    print "\n", "----- HERE IS THE LINEMERGE OUTPUT ------"
    wkt_print_geoms(newgeoms)
    
    ################
    # POLYGONIZE
    ################
    temp = geos.vector_GeometryP()
    for g in geoms:
        temp.push_back(g)
    plgnzr = geos.Polygonizer()
    plgnzr.add(temp)
    polys = plgnzr.getPolygons()
    newgeoms = []
    for i in range(polys.size()):
        newgeoms.append(polys[i])
    del polys

    print "\n", "----- HERE IS POLYGONIZE OUTPUT ------"
    wkt_print_geoms(newgeoms)

print "GEOS", geos.geosversion(), "ported from JTS", geos.jtsport()
do_all()