DM_numpy_spatial_query.py

## @package python.demo.DM_numpy_spatial_query
# "
# example script showing how to load an OPALS datamanager 
# and perform a spatial query returned as numpy dictionary
# 
 
from __future__ import print_function  # print function syntax as in python 3
import sys
 
from opals import pyDM
 
def DM_spatial_query(filename, renderPoints = False):
  # pyDM.Datamanager.load parameter: filename(string), readOnly(bool) threadSafety(bool)
  dm = pyDM.Datamanager.load(filename, True, False)
  if not dm:
    print("Unable to open ODM '" + filename + "'")
    sys.exit(1)
 
  # output some pyDM statistic
  limit = dm.getLimit()
  print("ODM contains", dm.sizePoint(), "points")
  print("2D-limit (%.3f," % limit.xmin, "%.3f) -" % limit.ymin, "(%.3f," % limit.xmax, "%.3f)" % limit.ymax)
 
  # build layout for attributes of interest (subset of odm attributes)
  lf = pyDM.AddInfoLayoutFactory()
  type, inDM = lf.addColumn(dm, "Id", True);             assert inDM == True
  type, inDM = lf.addColumn(dm, "GPSTime", True);        assert inDM == True
  type, inDM = lf.addColumn(dm, "Amplitude", True);      assert inDM == True
  type, inDM = lf.addColumn(dm, "Classification", True); assert inDM == True
  layout = lf.getLayout()
 
  # Create query object
  # polygon
  pf = pyDM.PolygonFactory()
  pf.addPoint(529600, 5338600)
  pf.addPoint(529650, 5338620)
  pf.addPoint(529650, 5338650)
  pf.addPoint(529630, 5338650)
  pf.closePart()
  queryPoly = pf.getPolygon()
  queryWin = pyDM.Window(529600, 5338600, 529650, 5338650) # window
 
 
  # creaate optional filter. select only ground points
  filter = pyDM.Filter("class[ground]")
 
  # [query points as dict of numpy arrays]
  print("\nPerform spatial query")
  result = pyDM.NumpyConverter.searchPoint(dm, queryPoly, layout, withCoordinates=True) # polygon query without filter
  #result = pyDM.NumpyConverter.searchPoint(dm, queryWin, layout, withCoordinates=True, filter=filter) # window query with filter
 
  print("result=", result)  # print returned object
  # [query points as dict of numpy arrays]
 
  # number of points found
  ptsQueried = result[next(iter(result))].size
  print(f"{ptsQueried} points queried")
 
  if renderPoints:
    # render points if requested
    import matplotlib.pyplot as plt
    fig = plt.figure()
    ax = fig.add_subplot(projection='3d')
    ax.scatter(result["x"], result["y"], result["z"], c=result["Amplitude"], marker='.', cmap=plt.cm.Greys.reversed())
 
    plt.axis('equal')
    plt.show()
 
 
if len(sys.argv) == 1:
    print("ODM parameter missing")
    sys.exit(-1)
 
renderPts = False
if len(sys.argv) > 2:
    renderPts = sys.argv[2] == '1' or sys.argv[2].lower() == 'true'
 
DM_spatial_query(sys.argv[1], renderPts)