[GIS] How to Convert 3D coordinates to 2D coordinates for R-Tree

coordinate systempythonshapelyspatial statistics

I am trying to implement a bounding box via a spatial index and get an error stating that I must either provide coordinates as either (minx, miny, maxx, maxy) or (x,y) for 2D indexes.

I've tried using the R-Tree example from this thread: More Efficient Spatial join in Python without QGIS, ArcGIS, PostGIS, etc

The first time I try to do a spatial index, I have the following value for a point ("Z" is the current value of "i" before my code breaks):

POINT Z (11.52951677300007 0.7729100360000416 -50000)

How can I gracefully get rid of the Z-coordinate to pass on as a 2D coordinate in my "for j" statement? Code snippet is below:

    from rtree import index
    idx = index.Index() #Create an R-Tree index and store the features in it (bounding box)
    for pos, poly in enumerate(polygons):
        idx.insert(pos, shape(poly['geometry']).bounds)

        for i,pt in enumerate(points):  #iterate through the points
            point = shape(pt['geometry'])
            for j in idx.intersection(point.coords[0]):  #iterate through spatial index
                if point.within(shape(multi[j]['geometry'])):
                    discoveredBuildings.append(point['id'])

Best Answer

It is a pure Python problem

With shapely, your 3D point is represented by

pt3D = Point(11.52951677300007,0.7729100360000416, -50000)
list(pt3D.coords)
[(11.52951677300007, 0.7729100360000416, -50000.0)]

Using slicing for example (there are others solutions)

pt2D = Point(list(pt3D.coords)[0][:2])
print pt2D
POINT (11.52951677300007 0.7729100360000416)

Or

 Pt2D = Point(pt3D.x, pt3D.y)
 print Pt2D
 POINT (11.52951677300007 0.7729100360000416)

Or ...