import numpy as np
pc = np.loadtxt('https://raw.githubusercontent.com/OSGeoLabBp/tutorials/master/english/data_processing/lessons/code/lidar.txt', delimiter=',') 
# load point cloud from text file

"""
Explanation: <a href="https://colab.research.google.com/github/OSGeoLabBp/tutorials/blob/master/english/data_processing/lessons/point_cloud_section.ipynb" target="_parent"><img src="https://colab.research.google.com/assets/colab-badge.svg" alt="Open In Colab"/></a>
Point cloud section
Getting either a horizontal or a vertical section of a point cloud is quite a common job in Geoinformatics these days. There is a wide variety of software to complete this task. For the very first trials, the open-source CloudCompare would be an option, having both graphical and command line interfaces available.
In this tutorial, some basic python scripts are presented in order to be able to set up your own preferences. 
First, load the point cloud into a numpy array. The file contains UTM easting and northing coordinates as well as elevation above mean sea level. Fields are separated by comma and there is no header.
End of explanation
"""


pc.shape[0], np.min(pc[:,2]), np.max(pc[:,2])

"""
Explanation: Check out the number of points loaded and get the minimum and maximum of the elevations
End of explanation
"""


tol = 1 # tolerance
coo = 1100 # elevation of section plane
sec = pc[np.absolute(pc[:, 2] - coo) < tol]

"""
Explanation: Getting a horizontal section is nothing else than just filtering out those points whose height differs from the section plane within a given tolerance.
End of explanation
"""


sec.shape[0]

"""
Explanation: Check out the number of points in the section:
End of explanation
"""


import matplotlib.pyplot as plt
plt.plot(sec[:,0], sec[:,1], 'o')
plt.xlabel('east (m)')
plt.ylabel('north (m)')
plt.axis('equal')
plt.grid('on')
plt.title('section at level ' + str(coo) + ' m')

"""
Explanation: and plot the points in a top view:
End of explanation
"""


e1 = 549400; n1 = 5128900 #first point of the section line
e2 = 549200; n2 = 5129300 #second point of the section line

"""
Explanation: Another common job worth discussing in this tutorial is the need for the creation of a vertical section of the point cloud along a line. Without going into the details, this job will be completed using homogenous coordinates. The section line is defined by the coordinates of two points on the line, as follows:
End of explanation
"""


from math import hypot
def vplain(x1, y1, x2, y2):
  """ set up line equation
  vp[0] * x + vp[1] * y + vp[2] = 0

  x1, y1 - horizontal coordinates of the start point of the section
  x2, y2 - horizontal coordinates of the end point of the section
  returns a numpy array with coefficients of the vertical plane
  """

  vp = np.zeros((3,))
  vp[0] = y1 - y2
  vp[1] = x2 - x1
  vp[2] = x1 * y2 - x2 * y1
  vp = vp / hypot(vp[0], vp[1]) # normalize
  return vp

vplain(e1, n1, e2, n2)

"""
Explanation: At first, we need to define a function to set up the line equation. In case a point is exactly on the line, the point coordinates fulfill the equations:
End of explanation
"""


def section(pc, x1, y1, x2, y2, tol):
  """ Select point from a point cloud near to a line

  pc - point cloud in a numpy array
  x1, y1 - horizontal coordinates of the start point of the section
  x2, y2 - horizontal coordinates of the end point of the section
  tol - tolerance distance from the section
  returns a numpy array with points near to the section
  """
  pc1 = pc.copy()
  pc1[:, 2] = 1 # change to homogenous coordinates
  vp = vplain(x1, y1, x2, y2) # equation of vertical plain
  sec = pc[np.abs(np.dot(pc1, vp)) < tol] # select points close to the section

  return sec

"""
Explanation: At second, we need to filter out those points whose distance from the section line is within the tolerance. The distance is computed by multiplying the point cloud and the array with the coefficients of the line. Strange as it may sound, but now the distance has a sign that informs whether a specific point is on the left or on the right side of the line. For filtering, we use the absolute value of the distances.
End of explanation
"""


sec = section(pc, e1, n1, e2, n2, tol)
sec.shape[0]

"""
Explanation: Check out the number of points in the section:
End of explanation
"""


plt.plot(sec[:,0], sec[:,1], 'o')
plt.xlabel('east (m)')
plt.ylabel('north (m)')
plt.axis('equal')
plt.grid('on')

"""
Explanation: Plot the points in a top view:
End of explanation
"""


def tr(e1, n1, e2, n2):
  """ set up transformation matrix for homogenous coordinates

  Parameters:
  e1, n1 - start point of the section line
  e2, n2 - end point of the section section line
  returns the transformation matrix
  """
  de = e2 - e1
  dn = n2 - n1

  d = hypot(de, dn) # distance
  r = de / d # sin
  m = dn / d # cos
  return np.dot(np.array([[1, 0, 0], [0, 1, 0], [-e1, -n1, 1]]),
                np.array([[m, r, 0], [-r, m, 0], [0, 0, 1]]))
  

"""
Explanation: In order to plot a section view, we need to transform section points into a system defined by the section line. To put it another way, the transformation matrix contains a shift to the first point of the section line and a rotation by the whole circle bearing of the line.
End of explanation
"""


# transformation matrix
trm = tr(e1, n1, e2, n2)
# make a copy of section points for homogenous transformation
pc1 = sec.copy()
pc1[:, 2] = 1
pc1 = np.dot(pc1, trm) # rotate points into the section plain
pc1[:, 2] = sec[:, 2] # copy back elevations to transformed points

"""
Explanation: Then apply the transformation to the section points:
End of explanation
"""


plt.plot(pc1[:,1], pc1[:,2], 'o')
plt.xlabel('chainage (m)')
plt.ylabel('elevation (m)')
plt.axis('equal')
plt.grid('on')

"""
Explanation: Plot the points in a section view:
End of explanation
"""


