<jupyter_start><jupyter_text>What is the proposed task:- ingest some liDAR points into a HDF file- ingest the aircraft trajectory into the file- anything else...and then **extract** data from the HDF file at different rates using a spatial 'query' What do the data look like?ASCII point clouds with the following attributes, currently used as a full M x N array:- time (GPS second of week, float)- X coordinate (UTM metres, float)- Y coordinate (UTM metres, float)- Z coordinate (ellipsoidal height, float)- return intensity (unscaled, float)- scan angle (degrees, float)Everything above is easily stored in the binary .LAS format (or .LAZ). It is kept in ASCII because the following additional data have no slots in .LAS:- X uncertainty (aircraft reference frame, m, float)- Y uncertainty (aircraft reference frame, m, float)- Z uncertainty (aircraft reference frame, m, float)- 3D uncertainty (metres, float)...and optionally (depending on the use case):- aircraft trajectory height (to ITRF08, metres)- aicraft position uncertainty X (metres, relative to aircraft position)- aircraft and sensor attributes...and derived data:- sea ice elevations (m, float)- estimted snow depths (m, float)- estimated snow depth uncertainty (m, float)- estimated ice thickness (m, float)- estimated ice thickness uncertainty (m, float)So, you can see how quickly .LAS loses it's value. ASCII point clouds are conceptually simple, but very big - and not well suited to use in a HPC context. Too much storage overhead, and you have to read the entire file in order to extract a subset. Six million points gets to around 50MB, it's pretty inefficient.So lets look at about 6 million points...  Let's look at a small set of pointsScenario:  A small set of LiDAR and 3D photogrammetry points collected adjacent to a ship (RSV Aurora Australia) parked in sea ice.Source: AAD's APPLS system (http://seaice.acecrc.org.au/crrs/)https://data.aad.gov.au/aadc/metadata/metadata.cfm?entry_id=SIPEX_II_RAPPLShttps://data.aad.gov.au/aadc/metadata/metadata.cfm?entry_id=SIPEX_LiDAR_sea_icePretty cover photo:<jupyter_code>import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
%matplotlib inline
#import plot_lidar
from datetime import datetime<jupyter_output><empty_output><jupyter_text>import a LiDAR swath<jupyter_code>swath = np.genfromtxt('../../PhD/python-phd/swaths/is6_f11_pass1_aa_nr2_522816_523019_c.xyz')
import pandas as pd

columns = ['time', 'X', 'Y', 'Z', 'I','A', 'x_u', 'y_u', 'z_u', '3D_u']
swath = pd.DataFrame(swath, columns=columns)
swath[1:5]<jupyter_output><empty_output><jupyter_text>Now load up the aircraft trajectory<jupyter_code>air_traj = np.genfromtxt('../../PhD/is6_f11/trajectory/is6_f11_pass1_local_ice_rot.3dp')
columns = ['time', 'X', 'Y', 'Z', 'R', 'P', 'H', 'x_u', 'y_u', 'z_u', 'r_u', 'p_u', 'h_u']
air_traj = pd.DataFrame(air_traj, columns=columns)
air_traj[1:5]<jupyter_output><empty_output><jupyter_text>take a quick look at the data<jupyter_code>fig = plt.figure(figsize = ([30/2.54, 6/2.54]))
ax0 = fig.add_subplot(111)      
a0 = ax0.scatter(swath['Y'],  swath['X'], c=swath['Z'] - np.min(swath['Z']), cmap = 'gist_earth',
                vmin=0, vmax=10, edgecolors=None,lw=0, s=0.6)
a1 = ax0.scatter(air_traj['Y'], air_traj['X'], c=air_traj['Z'], cmap = 'Reds',
                lw=0, s=1)
plt.tight_layout()<jupyter_output><empty_output><jupyter_text>Making a HDF file out of those points<jupyter_code>import h5py
#create a file instance, with the intention to write it out
lidar_test = h5py.File('lidar_test.hdf5', 'w')
swath_data = lidar_test.create_group('swath_data')
swath_data.create_dataset('GPS_SOW', data=swath['time'])
#some data
swath_data.create_dataset('UTM_X', data=swath['X'])
swath_data.create_dataset('UTM_Y', data=swath['Y'])
swath_data.create_dataset('Z', data=swath['Z'])
swath_data.create_dataset('INTENS', data=swath['I'])
swath_data.create_dataset('ANGLE', data=swath['A'])
swath_data.create_dataset('X_UNCERT', data=swath['x_u'])
swath_data.create_dataset('Y_UNCERT', data=swath['y_u'])
swath_data.create_dataset('Z_UNCERT', data=swath['z_u'])
swath_data.create_dataset('3D_UNCERT', data=swath['3D_u'])
#some attributes
lidar_test.attrs['file_name'] = 'lidar_test.hdf5'

lidar_test.attrs['codebase'] = 'https://github.com/adamsteer/matlab_LIDAR'<jupyter_output><empty_output><jupyter_text>That's some swath data, now some trajectory data at a different sampling rate<jupyter_code>traj_data = lidar_test.create_group('traj_data')
#some attributes
traj_data.attrs['flight'] = 11
traj_data.attrs['pass'] = 1
traj_data.attrs['source'] = 'RAPPLS flight 11, SIPEX-II 2012'
#some data
traj_data.create_dataset('pos_x', data = air_traj['X'])
traj_data.create_dataset('pos_y', data =  air_traj['Y'])
traj_data.create_dataset('pos_z', data =  air_traj['Z'])<jupyter_output><empty_output><jupyter_text>close and write the file out<jupyter_code>lidar_test.close()<jupyter_output><empty_output><jupyter_text>OK, that's an arbitrary HDF file builtThe generated file is substantially smaller than the combined sources - 158 MB from 193, with no attention paid to optimisation.The .LAZ version of the input text file here is 66 MB. More compact, but we can't query it directly - and we have to fake fields! Everything *in* the swath dataset can be stored, but we need to pretend uncertainties are RGB, so if person X comes along and doesn't read the metadata well, they get crazy colours, call us up and complain. Or we need to use .LAZ extra bits, and deal with awkward ways of describing things.It's also probably a terrible HDF, with no respect to CF compliance at all. That's to come :)  And now we add some 3D photogrammetry at about 80 points/m^2:<jupyter_code>photo = np.genfromtxt('/Users/adam/Documents/PhD/is6_f11/photoscan/is6_f11_photoscan_Cloud.txt',skip_header=1)
columns = ['X', 'Y', 'Z', 'R', 'G', 'B']
photo = pd.DataFrame(photo[:,0:6], columns=columns)
#create a file instance, with the intention to write it out
lidar_test = h5py.File('lidar_test.hdf5', 'r+')
photo_data = lidar_test.create_group('3d_photo')
photo_data.create_dataset('UTM_X', data=photo['X'])
photo_data.create_dataset('UTM_Y', data=photo['Y'])
photo_data.create_dataset('Z', data=photo['Z'])
photo_data.create_dataset('R', data=photo['R'])
photo_data.create_dataset('G', data=photo['G'])
photo_data.create_dataset('B', data=photo['B'])
#del lidar_test['3d_photo']
lidar_test.close()<jupyter_output><empty_output><jupyter_text>Storage is a bit less efficient here.- ASCII cloud: **2.1 GB**- .LAZ format with same data: **215 MB**- HDF file containing LiDAR, trajectory, 3D photo cloud: **1.33 GB**So, there's probably a case for keeping super dense clouds in different files (along with all their ancillary data). Note that .LAZ is able to store all the data used for the super dense cloud here. But - how do we query it efficiently?Also, this is just a demonstration, so we push on!  now, lets look at the HDF file... and get stuff<jupyter_code>from netCDF4 import Dataset
thedata = Dataset('lidar_test.hdf5', 'r')
thedata<jupyter_output><empty_output><jupyter_text>There are the two groups - swath_data and traj_data<jupyter_code>swath = thedata['swath_data']
swath
utm_xy = np.column_stack((swath['UTM_X'],swath['UTM_Y']))
idx = np.where((utm_xy[:,0] > -100) & (utm_xy[:,0] < 200) & (utm_xy[:,1] > -100) & (utm_xy[:,1] < 200)  )
chunk_z = swath['Z'][idx]
chunk_z.size
max(chunk_z)
chunk_x = swath['UTM_X'][idx]
chunk_x.size
chunk_y = swath['UTM_Y'][idx]
chunk_y.size
chunk_uncert = swath['Z_UNCERT'][idx]
chunk_uncert.size
plt.scatter(chunk_x, chunk_y, c=chunk_z, lw=0, s=2)<jupyter_output><empty_output><jupyter_text>That gave us a small chunk of LIDAR points, without loading the whole point dataset. Neat!...but being continually dissatisfied, we want more! Lets get just the corresponding trajectory:<jupyter_code>traj = thedata['traj_data']
traj<jupyter_output><empty_output><jupyter_text>Because there's essentiually no X extent for flight data, only the Y coordinate of the flight data are needed...<jupyter_code>pos_y = traj['pos_y']
idx = np.where((pos_y[:] > -100.) & (pos_y[:] < 200.))
cpos_x = traj['pos_x'][idx]
cpos_y = traj['pos_y'][idx]
cpos_z = traj['pos_z'][idx]<jupyter_output><empty_output><jupyter_text>Now plot the flight line and LiDAR together<jupyter_code>plt.scatter(chunk_x, chunk_y, c=chunk_z, lw=0, s=3, cmap='gist_earth')
plt.scatter(cpos_x, cpos_y, c=cpos_z, lw=0, s=5, cmap='Oranges')<jupyter_output><empty_output><jupyter_text>...and prove that we are looking at a trajectory *and* some LiDAR<jupyter_code>from mpl_toolkits.mplot3d import Axes3D
#set up a plot
plt_az=310
plt_elev = 40.
plt_s = 3
cb_fmt = '%.1f'

cmap1 = plt.get_cmap('gist_earth', 10)

#make a plot
fig = plt.figure()
fig.set_size_inches(35/2.51, 20/2.51)
ax0 = fig.add_subplot(111, projection='3d')
a0 = ax0.scatter(chunk_x, chunk_y, (chunk_z-min(chunk_z))*2,
                 c=np.ndarray.tolist((chunk_z-min(chunk_z))*2),\
                cmap=cmap1,lw=0, vmin = -0.5, vmax = 5, s=plt_s)
ax0.scatter(cpos_x, cpos_y, cpos_z, c=np.ndarray.tolist(cpos_z),\
                cmap='hot', lw=0, vmin = 250, vmax = 265, s=10)
ax0.view_init(elev=plt_elev, azim=plt_az)
plt.tight_layout()<jupyter_output><empty_output><jupyter_text>plot coloured by point uncertainty<jupyter_code>#set up a plot
plt_az=310
plt_elev = 40.
plt_s = 3
cb_fmt = '%.1f'

cmap1 = plt.get_cmap('gist_earth', 30)

#make a plot
fig = plt.figure()
fig.set_size_inches(35/2.51, 20/2.51)
ax0 = fig.add_subplot(111, projection='3d')
a0 = ax0.scatter(chunk_x, chunk_y, (chunk_z-min(chunk_z))*2,
                 c=np.ndarray.tolist(chunk_uncert),\
                 cmap=cmap1, lw=0, vmin = 0, vmax = 0.2, s=plt_s)
ax0.scatter(cpos_x, cpos_y, cpos_z, c=np.ndarray.tolist(cpos_z),\
                cmap='hot', lw=0, vmin = 250, vmax = 265, s=10)
ax0.view_init(elev=plt_elev, azim=plt_az)
plt.tight_layout()
plt.savefig('thefig.png')<jupyter_output><empty_output><jupyter_text>now pull in the photogrammetry cloudThis gets a little messy, since it appears we still need to grab X and Y dimensions - so still 20 x 10^6 x 2 points. Better than 20 x 10^6 x 6, but I wonder if I'm missing something about indexing.<jupyter_code>photo = thedata['3d_photo']
photo
photo_xy = np.column_stack((photo['UTM_X'],photo['UTM_Y']))
idx_p = np.where((photo_xy[:,0] > 0) & (photo_xy[:,0] < 100) & (photo_xy[:,1] > 0) & (photo_xy[:,1] < 100)  )
plt.scatter(photo['UTM_X'][idx_p], photo['UTM_Y'][idx_p], c = photo['Z'][idx_p],\
                cmap='hot',vmin=-1, vmax=1, lw=0, s=plt_s)
p_x = photo['UTM_X'][idx_p]
p_y = photo['UTM_Y'][idx_p]
p_z = photo['Z'][idx_p]
plt_az=310
plt_elev = 70.
plt_s = 2

#make a plot
fig = plt.figure()
fig.set_size_inches(25/2.51, 10/2.51)
ax0 = fig.add_subplot(111, projection='3d')

#LiDAR points
ax0.scatter(chunk_x, chunk_y, chunk_z-50, \
            c=np.ndarray.tolist(chunk_z),\
            cmap=cmap1, vmin=-30, vmax=2, lw=0, s=plt_s)

#3D photogrammetry pointd
ax0.scatter(p_x, p_y, p_z, 
            c=np.ndarray.tolist(p_z),\
            cmap='hot', vmin=-1, vmax=1, lw=0, s=5)

#aicraft trajectory
ax0.scatter(cpos_x, cpos_y, cpos_z, c=np.ndarray.tolist(cpos_z),\
                cmap='hot', lw=0, vmin = 250, vmax = 265, s=10)


ax0.view_init(elev=plt_elev, azim=plt_az)
plt.tight_layout()
plt.savefig('with_photo.png')<jupyter_output><empty_output><jupyter_text>This is kind of a clunky plot - but you get the idea (I hope). LiDAR is in blues, the 100 x 100 photogrammetry patch in orange, trajectory in orange. Different data sources, different resolutions, extracted using pretty much the same set of queries.<jupyter_code>print('LiDAR points: {0}\nphotogrammetry points: {1}\ntrajectory points: {2}'.
      format(len(chunk_x), len(p_x), len(cpos_x) ))<jupyter_output>LiDAR points: 43264
photogrammetry points: 817370
trajectory points: 1448
{'fraction_non_alphanumeric': 0.13896136688830907, 'fraction_numerical': 0.027775439777796482, 'mean_word_length': 5.539350577875619, 'pattern_counts': {'":': 0, '<': 73, '<?xml version=': 0, '>': 73, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '588743', 'license': 'Apache-2.0', 'n_tokens_mistral': 4736, 'n_tokens_neox': 4476, 'n_words': 1164, 'path': '.ipynb_checkpoints/Point cloud to HDF-checkpoint.ipynb', 'repo_name': 'adamsteer/nci-notebooks', 'size': '11881'}