<jupyter_start><jupyter_text>Site 1: <NAME>, September 24th 2019 This notebook is used to derive results presented in Section 3.2 of the paper:**Digitalization of scanning lidar measurement campaign planning**The results are generated using the Python library *campaign-planning-tool* v0.1.3.To be able to run this notebook user are requested to install *campaign-planning-tool*. Use the following link for installation instructions: https://github.com/niva83/campaign-planning-tool<jupyter_code>from campaign_planning_tool import CPT
import numpy as np
from random import shuffle
import matplotlib.pyplot as plt<jupyter_output><empty_output><jupyter_text>Let us start by creating a CPT object and configuring storage and landcover data paths:<jupyter_code>layout = CPT()
layout.set_path('../data/output/', path_type = 'output')
layout.set_path('../data/input/landcover/glc_V18_5.tif', path_type = 'landcover')<jupyter_output>Path /Users/niva/Desktop/supplementary_material set for storing CPT outputs
Path /Volumes/Secondary_Drive/work/projects/campaign-planning-tool/examples/../data/input/landcover/glc_V18_5.tif set for landcover data<jupyter_text>Before we proceed adding the measurement points it important to know that the CPT methods perform calculation in UTM coordinate system (i.e., positions provided as triplets Easting, Northing and Height). Therefore it is required to set a proper UTM zone to the class, which in our case is 30U. The digits in the UTM zone represent so-called latitudinal zone, while a character represent longitudinal zone. Both, digits and character are required to be provided to the CPT class:<jupyter_code>layout.set_utm_zone('30U')<jupyter_output>Correct latitudinal zone!
Correct longitudinal zone!
UTM zone set<jupyter_text>Now that we have set a proper UTM zone we can add our measurement point to the CPT class.In this example we use turbine hub positions as measurement points:<jupyter_code>points = np.array([
[534989,6184410,47],
[535205,6184536,47],
[535523,6184458,47],
[535837,6184425,47],
[536008,6184602,47],
[536288,6184670,47],
[535150,6184166,47],
[535468,6184187,47],
[535119,6183722,47],
[535382,6183935,47],
[535646,6183803,47],
[535865,6183938,47],
[536172,6184027,47],
[536404,6184103,47],
[535465,6183116,47],
[535432,6183632,47],
[535707,6183551,47],
[535958,6183688,47],
[536209,6183662,47],
[535636,6183310,47],
[535661,6182950,47],
[535440,6182744,47]
])<jupyter_output><empty_output><jupyter_text>Since CPT methods work with the height provided as the height above sea level, thus we need to add terrain height to the turbine hub height. The CPT class has a method to fetch the terrain height from SRTM database:<jupyter_code>points[:, 2] = points[:, 2] + layout.get_elevation('30U', points)<jupyter_output>Correct latitudinal zone!
Correct longitudinal zone!<jupyter_text>Now lets add the measurement points to the CPT class, specifically to the measurement point dictionary as *initial* measurements:<jupyter_code>layout.add_measurement_instances('initial', points)
layout.measurements_dictionary['initial']<jupyter_output><empty_output><jupyter_text>Before we start generating GIS layers which will help us navigate positioning of two lidars let us change the following attributes of the CPT class:- *maximum allowed elevation angle*- *expected average range of lidards*<jupyter_code>layout.MAX_ELEVATION_ANGLE = 5 # in degrees
layout.AVERAGE_RANGE = 3000 # in meters<jupyter_output><empty_output><jupyter_text>Now lets generate the GIS layer for placing of the first lidar. The method that we need to call is:```layout.generate_first_lidar_placement_layer(points_id)```where point_id is a string indicating which measurement points should be used in this process. We will set *points_id* to '*initial*' since we have previously added measurement points considering this id:<jupyter_code>layout.generate_first_lidar_placement_layer('initial')
layout.plot_layer('first_lidar_placement', title = 'First lidar placement map')<jupyter_output><empty_output><jupyter_text>The plot shows you areas where if you place lidar you will be able to reach a certain number of measurement points considering a number of constriants such as:* Range* Unobstracted line-of-sight* Elevation angle* Restriction zones Currently the CPT library don't provide an interactive way of placing lidar or optimization routine to do this for you. Therefore, for time being you will have to this manually.You have two options of performing this task. You can export the above layer to KML by using following built-in method:```layout.lidar_position_suggestion(filename, **kwargs)```You can use another built-in method to get a list of potential positions:```layout.lidar_position_suggestion(layer_id, teshold)```We will demonstrate both methods:<jupyter_code>layout.export_kml('first_lidar', layer_ids = ['first_lidar_placement'])
first_lidar_positions = layout.lidar_position_suggestion('first_lidar_placement', 22)<jupyter_output>KML exported successful!<jupyter_text>In the first method we are providing the KML file name followed with *layer_id* suppiled as a list.In the second method we provide *layer_id* string followed with a *treshold* provided as an integer.  The treshold represents a minimum number of measurement points the suggested lidar position should be able to provide considering all the constraints. In our case, we want only locations which will allow us to reach all 22 measurement points. Nevertheless, in this example, using Google Earth we have identified position for the first lidar and add it to the lidar dictionary:<jupyter_code>layout.add_lidar_instance('ws1', np.array([536062, 6182292, 350]), layer_id = 'first_lidar_placement')
layout.plot_layer('first_lidar_placement', lidar_ids = ['ws1'])<jupyter_output>Lidar 'ws1' added to the lidar dictionary, which now contains 1 lidar instance(s).
Updating lidar instance 'ws1' considering GIS layer 'first_lidar_placement'.<jupyter_text>What happened in the previous two lines of code? First we call method```layout.add_lidar_instance(lidar_id, position, kwargs)```which adds lidar position to the lidar dictionary.In our case we provided this method with following parameters:* lidar_id ='*ws1*'* position, which was identified by inspecting GIS layer in Google Earth* layer_id = '*first_lidar_placement*'The last parameter points the method to extract information about the reachable points from a specific GIS layer, which in the case of the first lidar placement is *'first_lidar_placement'* and updates the lidar dictionary instance 'ws1' accordingly. You can access the lidar dictionary instance by calling:```layout.lidar_dictionary['ws1']```, and if you inspect specifically the following key:```layout.lidar_dictionary['ws1']['reachable_points']```you can see that the returned array contains 22 elements (same lenght as points) and all elements are equal to 1 since all measurement points are reachable. Otherwise, if a point is unreachable by the lidar in the dictionary this will be indicated with 0. Inspect other elements of the lidar dictionary. We will use them at the later stage.<jupyter_code>layout.lidar_dictionary['ws1']['reachable_points']<jupyter_output><empty_output><jupyter_text>Since we have selected the first lidar position, let's add the second lidar. Now we need to add one more contraint, which is the minimum intersecting angle between laser beams.By calling ```layout.MIN_INTERSECTING_ANGLE```you can see what is the preset value in degrees.For time being we will use the default value of 30 degrees and create so called *'additional lidar placement layer'* To do this we will call method ```layout.generate_additional_lidar_placement_layer(lidar_id)```where *lidar_id* will be set to *'ws1'*, since that is the id of our first lidar.<jupyter_code>layout.generate_additional_lidar_placement_layer('ws1')
layout.plot_layer('additional_lidar_placement', lidar_ids = ['ws1'])<jupyter_output>Updating lidar instance 'ws1' considering GIS layer 'first_lidar_placement'.<jupyter_text>In the above plot of the new GIS layer you can see that some areas which in the previous plot were indicated as a good areas for lidar installation are removed. The reason for this is that if the second lidar was placed at those areas the intersecting angle between the first and second lidar would be lower than 30 degrees.Again we can export this newly made layer to KML file or use the internal method to suggest us the best positions for the second lidar:<jupyter_code>layout.export_kml('second_lidar', layer_ids = ['additional_lidar_placement'])
second_lidar_positions = layout.lidar_position_suggestion('additional_lidar_placement', 22)
layout.add_lidar_instance('ws2', np.array([537262, 6183492, 304]),layer_id = 'additional_lidar_placement')  
layout.plot_layer('additional_lidar_placement', lidar_ids = ['ws1', 'ws2'])<jupyter_output>Lidar 'ws2' added to the lidar dictionary, which now contains 2 lidar instance(s).
Updating lidar instance 'ws2' considering GIS layer 'additional_lidar_placement'.<jupyter_text>At this point we have positions of our dual-Doppler system and we are ready to optimize trajectory, plot campaign design and export lidars configurations. We will use the following built-in methods for this:```layout.optimize_trajectory(lidar_ids, **kwargs)layout.plot_design(layer_id, lidar_ids, **kwargs)layout.export_measurement_scenario(layer_id)```Additionally we will export KML containing layers, lidar positions and trajectory:```layout.export_kml(filename, **kwargs)```<jupyter_code>layout.optimize_trajectory(['ws1', 'ws2'], sync = True, only_common_points = True)<jupyter_output>Lidar instances:['ws1', 'ws2'] will be updated with the optimized trajectory
Synchronizing trajectories for lidar instances:['ws1', 'ws2']<jupyter_text>In this method we set kwargs *sync* and *only_common_points* to True. The first kwarg will assure synchronized trajectories, while the second kwarg assure that trajectory will consider only measurement points reachable by both lidars. The optimize_trajectory method currently only generates step-stare trajectories. You can access lidar dictionary and the result of this method (i.e., access key 'motion_config').<jupyter_code>layout.lidar_dictionary['ws1']['motion_config']
layout.lidar_dictionary['ws2']['motion_config']
layout.lidar_dictionary['ws2']['trajectory']
print('Trajectory motion time [s]: ' 
      + str(round(np.sum(layout.lidar_dictionary['ws1']['motion_config'].values[:,-1])/1000,2)))
layout.plot_design('additional_lidar_placement', lidar_ids = ['ws1', 'ws2'])
layout.export_measurement_scenario(['ws1', 'ws2'])
layout.export_kml('campaign_desing', 
                  layer_ids = ['first_lidar_placement','additional_lidar_placement'],
                  lidar_ids = ['ws1', 'ws2'])<jupyter_output>Measurement scenario export successful!
KML exported successful!<jupyter_text>The last method will export following files:1. Motion program to drive scanner heads (PMC file)2. Range gate file to configure laser and FPGA (TXT file)3. YAML and XML files containing human and machine readable compilation of information from (1) and (2)Motion programs and range gate files are currently only applicable for [long-range WindScanners](https://www.mdpi.com/2072-4292/8/11/896).Now you are all equiped to make scanning lidar measurements!  Appendix: verification of the optimized trajectoryFor the purpose of verifying that the optimized trajectory is really the shortest in duration we would need to generate all possible trajectories by making permutation of trajectory points. In this example this is not feasible since the trajectory contains 22 measurement points, thus there is 22! = 10^21 possible trajectories. Instead, we will make generate randomly 10^6 trajectories and calculate their timing.<jupyter_code>trajectory_init = layout.trajectory.values[:,1:].tolist()
lidar_ids= ['ws1', 'ws2']
i = 10**6
no_it = i
j = 0
sync_time_list = []
while(i>0):
    i = i - 1
    j = j + 1
    
    shuffle(trajectory_init)
    trajectory = np.asarray(trajectory_init)
    
    sync_time = []
    for lidar in lidar_ids:

        motion_table = layout.generate_trajectory(
                layout.lidar_dictionary[lidar]['position'], 
                trajectory)

        timing = motion_table.loc[:, 'Move time [ms]'].values
        sync_time = sync_time + [timing]
        
    # Synchronize the trajectory time for two lidars
    # by taking the maximum move time between two of them
    # for each trajectory point:        
    sync_time = np.sum(np.max(np.asarray(sync_time).T, axis = 1))
    sync_time_list = sync_time_list + [sync_time]   
    print("Progress {:2.1%}".format(j / no_it), end="\r")
sync_time_list = np.asarray(sync_time_list)
print('Average time: \t' + str(round(np.mean(sync_time_list)/1000,2)) +' s')
print('Max time: \t' + str(round(np.max(sync_time_list)/1000,2)) +' s')
print('Min time: \t' + str(round(np.min(sync_time_list)/1000,2)) +' s')
print('Std: \t\t' + str(round(np.std(sync_time_list)/1000,2)) +' s')


fig, ax = plt.subplots(sharey = True, 
                        figsize=(500/layout.MY_DPI, 
                                500/layout.MY_DPI), 
                        dpi=layout.MY_DPI)
n, bins, patches = plt.hist(sync_time_list/1000, 100, facecolor='black', alpha=1)
plt.xlabel('Trajectory time (motion only) [s]', fontsize = layout.FONT_SIZE)
plt.ylabel('Occurrence []', fontsize = layout.FONT_SIZE)

plt.show()<jupyter_output>Average time: 	22.09 s
Max time: 	25.76 s
Min time: 	17.15 s
Std: 		1.03 s
{'fraction_non_alphanumeric': 0.08993576017130621, 'fraction_numerical': 0.04009451377095178, 'mean_word_length': 5.892620865139949, 'pattern_counts': {'":': 0, '<': 56, '<?xml version=': 0, '>': 57, 'https://': 2, 'lorem ipsum': 0, 'www.': 1, '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': '16', 'dir': 'jupyter-structured-clean-dedup', 'id': '81726', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 4438, 'n_tokens_neox': 3989, 'n_words': 1408, 'path': 'site_1_scotland.ipynb', 'repo_name': 'niva83/campaign-planning-tool-examples', 'size': '13564'}