import numpy as np
import pandas as pd

#matplotlib.use('Agg')

import matplotlib
from matplotlib.collections import PolyCollection
import mpl_toolkits.mplot3d as a3
import matplotlib.colors as colors
from matplotlib import cm
import pylab as plt
from mpl_toolkits.mplot3d import Axes3D

import shapely.geometry as geometry
from descartes import PolygonPatch
from shapely.ops import cascaded_union, polygonize, unary_union, split
import math

%matplotlib inline
points = pd.read_csv("pufferfish1by1sample.xyz",
                       delim_whitespace=True,
                       names=["x","y","z","r","g","b"])
##### 1 meter by 1 meter sample at Pufferfish

print(f"1 meter sample at Pufferfish contains {len(points)} points")
print(points.head())

fig, ax = plt.subplots()
ax.plot(points['x'],points['y'],'.',color="slategrey")
ax.grid()
#fig.savefig("all.png")
plt.show()
1 meter sample at Pufferfish contains 496869 points
         x        y        z    r   g   b
0  10.6667 -12.4826 -18.9677   33   9   2
1  10.6669 -12.4816 -18.9680   41  21  13
2  10.6679 -12.4820 -18.9689   37  16  14
3  10.6683 -12.4811 -18.9689   61  43  41
4  10.6702 -12.4907 -18.9629  105  83  68

png

##### 50 cm sub-sample

samp = points.query("10.5 < x < 11.0")
samp = samp.query("-12.0 < y < -11.5")#.sample(10000)

print("\n")
print("50 cm sub-sample contains {} points.".format(len(samp)))
print("\n")
print(samp.head())
print("\n")

# fig = plt.figure()
# ax = fig.add_subplot(111, projection='3d')
# ax.scatter(samp['x'],samp['y'],samp['z'],marker='.',c="navy",alpha=0.5)
# ax.grid()
# plt.show()
50 cm sub-sample contains 112554 points.


             x        y        z    r    g   b
51610  10.5405 -11.9993 -18.9334  109  105  92
51613  10.5416 -11.9985 -18.9331  105  100  88
51614  10.5445 -11.9991 -18.9328  101   95  82
51615  10.5467 -11.9992 -18.9333  104   96  81
51616  10.5448 -11.9978 -18.9332  104   95  81

png

##### alternate solution technique
##### 1 cm slices

inc = 0.01 # 1 cm slices
lower = -12.0 # y minimum
upper = -11.5 # y maximum
AVF = []
shape_list = []

# Envelope size determined by entire 50 cm sample

lowerz = min(samp['z'])
upperz = max(samp['z'])
leftx = min(samp['x'])
rightx = max(samp['x'])
fronty = min(samp['y'])
backy = max(samp['y'])

envelope = geometry.Polygon([(leftx,lowerz),
                             (rightx,lowerz),
                             (rightx,upperz),
                             (leftx,upperz)])

for i in np.arange(lower,upper,inc):
    inc_slice = samp.query(f"{np.round(i,2)} < y <= {np.round(i+inc,2)}")

#     print(f"Slice for y = {np.round(lower,2)}, {np.round(lower+inc,2)} "+
#           f"contains {len(inc_slice)} points\n")

    minx = min(inc_slice['x'])
    maxx = max(inc_slice['x'])
    minz = min(inc_slice['z'])
    maxz = max(inc_slice['z'])

    # determine the left and right most points in slice
    left = (0,0)
    right = (0,0)
    for p in zip(inc_slice['x'],inc_slice['z']):
        if p[0] <= minx:
            left = (p[0],p[1])
        if p[0] >= maxx:
            right = (p[0],p[1])

    sliver_inc = 0.01 # 1 cm

##### Future Improvement
#     inc_pts = geometry.MultiPoint(list(zip(inc_slice['x'],inc_slice['y'],inc_slice['z'])))
#     envelope2 = inc_pts.envelope
#     print(envelope2)

    solid_area = [right,
                  (rightx,lowerz),
                  (leftx,lowerz),
                  left]

    for ix in np.arange(minx,maxx,sliver_inc):
        sliver = inc_slice.query(f"{ix} < x <= {ix+sliver_inc}")
        if len(sliver) == 0:
            new_point = (ix+sliver_inc/2.0,solid_area[-1][1])
#            print(f"{len(sliver)} points, {new_point}")
        else:
            new_point = (ix+sliver_inc/2.0,max(sliver['z']))
#            print(f"{len(sliver)} points, {new_point}")
        solid_area.append(new_point)

#    print(f"The area of the envelope is {envelope.area}")
    solid = geometry.Polygon(solid_area)
    shape_list.append(solid_area)
#    print(f"The area of the solid is {solid.area}")

    AVF.append(solid.area/envelope.area)
#    print(f"The solid area fraction for this slice is {AVF}")
print("done!")
done!
print(len(AVF))
print(np.mean(AVF))
50
0.4860697629701244
# plotting routine
fig = plt.figure()
ax = fig.add_subplot(111)
ax.add_patch(PolygonPatch(envelope,
                          fc='lightgrey',
                          ec='black',
                          alpha=0.5))

ax.scatter(inc_slice['x'],inc_slice['z'],marker='.',c="indigo")
ax.add_patch(PolygonPatch(solid,
                          fc='blue',
                          ec='black',
                          alpha=0.5))
ax.grid()

png

# plot of all slices
# y values on plot are exaggerated to provide spacing

facecolors = [cm.brg(x) for x in np.random.rand(50)]
zs = range(0,100,2)
fig = plt.figure(figsize=(11,8.5))
ax = fig.gca(projection='3d')
poly = PolyCollection(shape_list,
                      facecolors=facecolors,
                      edgecolor='k',
                      alpha=0.6)
#poly.set_color(colors.rgb2hex(sp.rand(3)))
ax.add_collection3d(poly,zs=zs,zdir='y')
ax.set_xlabel('X (m)')
ax.set_ylabel('False Y')
ax.set_zlabel('Z (m)')
ax.set_xlim(leftx, rightx)
ax.set_ylim(-1, 100)
ax.set_zlim(lowerz, upperz)

plt.show()

png