from bisect import bisect_left,bisect_right
from math import hypot,fabs
from p_ggen import iterby2
from p_gmath import thanSegSeguw
class ThanDTMlines:
"A set of lines which behaves as a Digital Terrain Model."
def __init__(self, dxmax=20.0, dext=50.0):
"Initialize DEM."
self.thanLines = []
self.thanDxmax = dxmax
self.thanDext = dext
def thanAddLine1(self, cp):
"Add a (3d) line to the DEM."
for ca, cb in iterby2(cp):
if ca[0] > cb[0]: ca, cb = cb, ca
dx = cb[0]-ca[0]
if dx <= self.thanDxmax:
self.thanLines.append((tuple(ca), tuple(cb)))
continue
n = int(dx/self.thanDxmax) + 1
ddx = self.thanDxmax/n
x = 0.0
cc = ca
for i in xrange(n-1):
x += ddx
cd = [ca1+(cb1-ca1)/dx*x for (ca1,cb1) in zip(ca, cb)]
self.thanLines.append((tuple(cc), tuple(cd)))
cc = cd
self.thanLines.append((tuple(cc), tuple(cb)))
def thanRecreate(self):
"Recreate DEM after lines additions."
if len(self.thanLines) < 1: return False, "DTM is empty!"
self.thanLines.sort()
return True, None
def __intersegZ(self, ca, cb):
"Compute intersections of segment with DEM lines; sort intersection from ca to cb."
ca = tuple(ca)
cb = tuple(cb)
# print "ca=", ca, "cb=", cb
rev = False
if ca[0] > cb[0]: ca, cb = cb, ca; rev=True
# print "ca=", ca, "cb=", cb, "rev=", rev
i = bisect_left(self.thanLines, ((ca[0]-self.thanDxmax,), (-1.0e30,)))
j = bisect_right(self.thanLines, ((cb[0],),(+1.0e30,)))
# print "i=", i, "j=", j, "len=", len(self.thanLines)
# for ii in xrange(i, j):
# print "thanlines", ii,":", self.thanLines[ii]
cint = []
for k in xrange(i, j):
c1, c2 = self.thanLines[k]
uw = thanSegSeguw(ca, cb, c1, c2)
if uw == None: continue
cd = [ca1+(cb1-ca1)*uw[1] for (ca1,cb1) in zip(c1, c2)]
if rev: cint.append((1-uw[0], cd))
else: cint.append((uw[0], cd))
cint.sort()
i = 0
j = 1
while j < len(cint):
if cint[j][0]-cint[i][0] < 0.001: # Delete double entries (with 0.1 percent anoxh)
del cint[j]
else:
i = j
j += 1
return cint
def thanPointZold(self, cp):
"Calculate the z coordinate of a point."
cp = list(cp) # bisect (called in __z1) needs cp as a list
ca = list(cp)
cb = list(cp)
ca[0] -= self.thanDext # |cb-ca| (1)
cb[0] += self.thanDext
d1, z1 = self.__z1(cp, ca, cb) # Try to find intersections along x direction
ca = list(cp)
cb = list(cp)
ca[1] -= self.thanDext # |cb-ca| must be the same as in (1), or the d1, d2..
cb[1] += self.thanDext # ..returned by __z1 will not be comparable
d2, z2 = self.__z1(cp, ca, cb) # Try to find intersections along y direction
if z1 == None: return z2 # Note that z2 may be None
if z2 == None: return z1
if d1 < d2: return z1
else: return z2
def thanPointZ(self, cp):
"Calculate the z coordinate of a point."
for p in (1, 2, 4):
cp = list(cp) # bisect (called in __z1) needs cp as a list
ca = list(cp)
cb = list(cp)
ca[0] -= self.thanDext*p # |cb-ca| (1)
cb[0] += self.thanDext*p
d1, z1 = self.__z1(cp, ca, cb) # Try to find intersections along x direction
# if z1 != None and fabs(z1) < 1.0: stop()
ca = list(cp)
cb = list(cp)
ca[1] -= self.thanDext*p # |cb-ca| must be the same as in (1), or the d1, d2..
cb[1] += self.thanDext*p # ..returned by __z1 will not be comparable
d2, z2 = self.__z1(cp, ca, cb) # Try to find intersections along y direction
# if z2 != None and fabs(z2) < 1.0: stop()
if z1 == None: z = z2 # Note that z2 may be None
elif z2 == None: z = z1
elif d1 < d2: z = z1
else: z = z2
if z != None: break
# if z != None and fabs(z) < 1.0: stop()
return z
def __z1(self, cp, ca, cb):
"Calclulate z with cp between ca and cb."
cint = self.__intersegZ(ca, cb)
# print "cint=", cint, "len=", len(cint)
if len(cint) < 1: return None, None
if len(cint) < 2:
d1 = fabs(cint[0][0]-0.5)
if d1 < 0.001: return d1, cint[0][1][2] # 1 intersection, but exactly on the isoline, with 0.1% tolerance
return None, None
if cint[0][0] > 0.5:
if cint[0][0] < 0.501: return cint[0][0]-0.5, cint[0][1][2] # Allow for 0.1% error
return None, None
if cint[-1][0] < 0.5:
if cint[-1][0] > 0.499: return 0.5-cint[-1][0], cint[-1][1][2] # Allow for 0.1% error
return None, None
j = bisect_left(cint, (0.5, cp))
# print "cint[j], j=", j
i = j-1
di = cint[i][0]
zi = cint[i][1][2]
dj = cint[j][0]
zj = cint[j][1][2]
z1 = zi+(zj-zi)/(dj-di)*(0.5-di)
d1 = min((0.5-di, dj-0.5))
# print "z1=", z1, "d1=", d1
return d1, z1
def thanLineZendpointstoo(self, cp):
"Calculate the z coordinates along the line cp."
cn = []
for ca, cb in iterby2(cp):
c = list(ca)
c[2] = self.thanPointZ(ca) # May be None
cn.append(c)
for d, c in self.__intersegZ(ca, cb):
if d < 0.001: continue # Avoid points near ca
if d > 0.999: break # Avoid points near cb
cn.append(c)
c = list(cb)
c[2] = self.thanPointZ(cb) # May be None
cn.append(c)
ni = interpolatez(cn)
return ni, cn
def thanLineZ(self, cp):
"Calculate the z coordinates along the line cp."
cn = []
for ca, cb in iterby2(cp):
c = list(ca)
c[2] = None # May be None
cn.append(c)
for d, c in self.__intersegZ(ca, cb):
if d < 0.001: continue # Avoid points near ca
if d > 0.999: break # Avoid points near cb
cn.append(c)
c = list(cb)
c[2] = self.thanPointZ(cb) # May be None
cn.append(c)
cn[0][2] = self.thanPointZ(cn[0])
ni = interpolatez(cn)
return ni, cn
def __interp2(cp):
"Interpolate the z coordinate to the points which have none."
d = [0.0]
j = len(cp) - 1
for k in xrange(1, j+1):
d.append(d[-1] + hypot(cp[k][1]-cp[k-1][1], cp[k][0]-cp[k-1][0]))
dij = d[-1]
zj = cp[j][2] # We access list in order, so we are fast
zi = cp[0][2]
fact = (zj - zi)/d[j]
for k in xrange(1, j):
cp[k][2] = zi + fact*d[k]
def interpolatez(cp):
"Interpolate the z coordinate to the points which have none; return number of interpolations."
n = len(cp)
for j in xrange(n):
if cp[j][2] != None: break
else:
return -1 # No z at all!
zj = cp[j][2]
for i in xrange(j): # Points from 0 to j-1 do not have valid z
cp[i][2] = zj
ni = j # Number of interpolated points
while True:
for i in xrange(j+1, n):
if cp[i][2] == None: break
else:
return ni # All points, from j to end, have valid z
for j in xrange(i+1, n):
if cp[j][2] != None: break
else:
break # No point, from j to end, has valid z
__interp2(cp[i-1:j+1])
ni += j-i-1
zj = cp[i-1][2]
for j in xrange(i, n):
cp[j][2] = zj
ni += 1
return ni
|