--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/safe.py	Sat May 29 11:07:34 2021 +0100
@@ -0,0 +1,121 @@
+#!/usr/bin/python3
+'''Find a safe landing space (i.e. in air, not above lava) near a cluster'''
+import sys, math
+
+from util import *
+from bisect import insort
+
+def usage():
+  print("""Usage:  safe.py [-d n] [filename]
+  n is maximum distance of landing to some point in the cluster, default is 10.
+  
+  Feed stdin with lines from the output of cluster.py, having produced
+  air.tsv with the same settings as the cluster input, and lava.tsv
+  with a lower-bound on Y of 0""")
+
+n=10.0
+if len(sys.argv)>1:
+  if sys.argv[1]=='-d':
+    sys.argv.pop()
+    n=float(sys.argv.pop(1))
+
+if len(sys.argv)>1:
+  infile=open(sys.argv.pop(1))
+else:
+  infile=sys.stdin
+
+class Block:
+  '''Shared by Air and Lava,
+  holds points and columns.
+  Columns are dictionaries indexed by (x,z) coords,
+  value is upward-ordered list of pairs, closed intervals of y'''
+  def __init__(self,filename):
+    with open(filename,'r') as file:
+      self.readHeaders(file)
+      self.readPoints(file)
+
+    cc=self.columns={}
+    tc={}
+    for p in self.points:
+      k=(p[0],p[2])
+      y=p[1]
+      try:
+        yy=tc[k]
+        insort(yy,y)
+      except KeyError:
+        tc[k]=[y]
+    # convert to lists of intervals
+    for k,yy in tc.items():
+      tyy=[yy]
+      ii=[]
+      while True:
+        clean=True
+        for j,yy in enumerate(tyy):
+          if len(yy)==1+(yy[-1]-yy[0]):
+            ii.append((yy[0],yy[-1]))
+          else:
+            clean=False
+            for i in range(len(yy)-1):
+              if yy[i]+1!=yy[i+1]:
+                ii+=(yy[0],yy[i])
+                tyy=[yy[i+1:]]+tyy[j+1:]
+            break
+        if clean:
+          break
+      self.columns[k]=set(ii) # so we can merge later
+
+  def readPoints(self,file):
+    file.readline()
+    self.points=[[float(i) for i in l.split('\t')[2].split(',')] for l in file]
+
+Block.readHeaders=readHeaders # from util
+
+class Air(Block):
+  pass
+
+A=Air('air.tsv')
+
+class Lava(Block):
+  pass
+
+L=Lava('lava.tsv')
+
+def d(p1,p2):
+  dx=p1[0]-p2[0]
+  dz=p1[1]-p2[1]
+  dy=p1[2]-p2[2]
+  return math.sqrt((dx*dx)+(dy*dy)+(dz*dz))
+
+def safety(p):
+  ka=[((a[0],a[2]),d(a,p)) for a in A.points if d(a,p)<=n]
+  return [(k,A.columns[k[0]],L.columns[k[0]]) for k in ka if k[0] in L.columns]
+
+readHeaders(sys.modules['__main__'],infile,False)
+for l in infile:
+  c=eval(l)
+  ss=[]
+  hits={}
+  misses=[]
+  for p in c:
+    ss=safety(p)
+    if ss:
+      for (k,a,l) in ss:
+        try:
+          (aa,ll)=hits[k]
+          aa.update(a)
+          ll.update(l)
+        except KeyError:
+          hits[k]=(a,l)
+    else:
+      misses.append(p)
+  print(c)
+  if hits:
+    print(' %s nearby landing columns'%len(hits))
+    for k in sorted(hits.keys(),key=lambda k:k[1]):
+      aa,ll=hits[k]
+      print("%s@%0.1f: %s %s"%(k[0],k[1],sorted(list(aa)),sorted(list(ll))))
+  else:
+    print(' No nearby landing zones')
+
+if infile!=sys.stdin:
+  infile.close()