add a rough visualiztion script that uses OGR to create a shapefile o…

…f the leaves of an index

scripts/visualize.py

@@ -0,0 +1,162 @@
+#!/usr/bin/env python
+
+from rtree import index
+
+import ogr
+
+
+def quick_create_layer_def( lyr, field_list):
+    # Each field is a tuple of (name, type, width, precision)
+    # Any of type, width and precision can be skipped.  Default type is string.
+
+    for field in field_list:
+        name = field[0]
+        if len(field) > 1:
+            type = field[1]
+        else:
+            type = ogr.OFTString
+
+        field_defn = ogr.FieldDefn( name, type )
+
+        if len(field) > 2:
+            field_defn.SetWidth( int(field[2]) )
+
+        if len(field) > 3:
+            field_defn.SetPrecision( int(field[3]) )
+
+        lyr.CreateField( field_defn )
+
+        field_defn.Destroy()
+
+
+import sys
+
+shape_drv = ogr.GetDriverByName('ESRI Shapefile')
+
+shapefile_name = sys.argv[1].split('.')[0]
+shape_ds = shape_drv.CreateDataSource( shapefile_name)
+leaf_block_lyr = shape_ds.CreateLayer( 'leaf', geom_type = ogr.wkbPolygon )
+point_block_lyr = shape_ds.CreateLayer( 'point', geom_type = ogr.wkbPolygon )
+point_lyr = shape_ds.CreateLayer( 'points', geom_type = ogr.wkbPoint )
+
+quick_create_layer_def( leaf_block_lyr,
+                                    [ 
+                                      ('BLK_ID', ogr.OFTInteger),
+                                      ('COUNT', ogr.OFTInteger),
+                                      ] )
+
+quick_create_layer_def( point_block_lyr,
+                                    [ 
+                                      ('BLK_ID', ogr.OFTInteger),
+                                      ('COUNT', ogr.OFTInteger),
+                                      ] )
+
+quick_create_layer_def( point_lyr,
+                                    [ 
+                                      ('ID', ogr.OFTInteger),
+                                      ('BLK_ID', ogr.OFTInteger),
+                                      ] )
+
+p = index.Property()
+p.filename = sys.argv[1]
+p.overwrite=False
+
+p.storage = index.RT_Disk
+idx = index.Index(sys.argv[1])
+
+
+leaves = idx.leaves()
+#leaves[0] == (0L, [2L, 92L, 51L, 55L, 26L], [-132.41727847799999, -96.717721818399994, -132.41727847799999, -96.717721818399994])
+
+from liblas import file
+f = file.File(sys.argv[1])
+
+def area(minx, miny, maxx, maxy):
+    width = abs(maxx - minx)
+    height = abs(maxy - miny)
+
+    return width*height
+
+def get_bounds(leaf_ids, lasfile, block_id):
+    # read the first point and set the bounds to that
+
+    p = lasfile.read(leaf_ids[0])
+    minx, maxx = p.x, p.x
+    miny, maxy = p.y, p.y
+
+    print len(leaf_ids)
+    print leaf_ids[0:10]
+
+    for p_id in leaf_ids:
+        p = lasfile.read(p_id)
+        minx = min(minx, p.x)
+        maxx = max(maxx, p.x)
+        miny = min(miny, p.y)
+        maxy = max(maxy, p.y)
+        feature = ogr.Feature( feature_def = point_lyr.GetLayerDefn() )
+        g = ogr.CreateGeometryFromWkt('POINT (%.8f %.8f)'%(p.x, p.y))
+        feature.SetGeometry(g)
+        feature.SetField('ID',p_id)
+        feature.SetField('BLK_ID',block_id)
+        result = point_lyr.CreateFeature(feature)
+
+    return (minx, miny, maxx, maxy)
+
+def make_poly(minx, miny, maxx, maxy):
+    wkt = 'POLYGON ((%.8f %.8f, %.8f %.8f, %.8f %.8f, %.8f %.8f, %.8f %.8f))' % (minx, miny, maxx, miny, maxx, maxy, minx, maxy, minx, miny)
+    shp = ogr.CreateGeometryFromWkt(wkt)
+    return shp
+
+def make_feature(lyr, geom, id, count):
+    feature = ogr.Feature( feature_def = lyr.GetLayerDefn() )
+    feature.SetGeometry(geom)
+    feature.SetField('BLK_ID',id)
+    feature.SetField('COUNT',count)
+    result = lyr.CreateFeature(feature)
+
+t = 0
+for leaf in leaves:
+    id = leaf[0]
+    ids = leaf[1]
+    count  = len(ids)
+    # import pdb;pdb.set_trace()
+
+    if len(leaf[2]) == 4:
+        minx, miny, maxx, maxy = leaf[2]
+    else:
+        minx, miny, maxx, maxy, minz, maxz = leaf[2]
+
+    if id == 186:
+        print leaf[2]
+
+
+
+    print leaf[2]
+    leaf = make_poly(minx, miny, maxx, maxy)
+    print 'leaf: ', minx, miny, maxx, maxy
+
+    pminx, pminy, pmaxx, pmaxy = get_bounds(ids, f, id)
+    point = make_poly(pminx, pminy, pmaxx, pmaxy)
+
+    print 'point:' ,pminx, pminy, pmaxx, pmaxy
+    print 'point bounds: ',point.GetArea(), area(pminx, pminy, pmaxx, pmaxy)
+    print 'leaf bounds: ', leaf.GetArea(), area(minx, miny, maxx, maxy)
+    print 'leaf - point: ', abs(point.GetArea() - leaf.GetArea())
+    print minx, miny, maxx, maxy
+  #  if shp2.GetArea() != shp.GetArea():
+  #      import pdb;pdb.set_trace()
+   # sys.exit(1)
+
+    make_feature(leaf_block_lyr, leaf, id, count)
+    make_feature(point_block_lyr, point, id, count)
+
+    t+=1
+   # if t ==2:
+#        break
+
+
+
+leaf_block_lyr.SyncToDisk()
+point_lyr.SyncToDisk()
+
+shape_ds.Destroy()