Data Processing Archive
This page is maintained as an archive of older methods. It was archived on 03/10/14. It contains references to the use of Mikrokopter flight controllers (ca. 2010) and the use of Data Processing files and scripts that have since been replaced by the tools available at code.ecosynth.org/ecosynthAerial
The Ecosynth data processing pipeline covers the steps necessary for preparing for a data collection mission, processing images into a 3D point cloud, then post-processing the point cloud data into useful data products like canopy height models (CHMs), digital terrain models (DTMs), maps of forest biomass density, and more...
Set up a GIS of the area of interest
Use existing high resolution imagery if available and a good DEM. If none are available then Google Earth high-resolution imagery is also useful and you can georeference a Google Earth image into GIS. Set up the following shapefiles:
- interest points - potential flying locations
- waypoints - waypoints to send to the hexa
- flight line - to build a uniform flight path
- area polygons - to identify maxi flying area, 250m radius, etc.
Design Flight Plan in GIS
- Currently Hexakopter can only go in a 250m radius from the point of take off so that the center of the AOI is usually the best place to launch from
- Tests have shown that an approx. 3.5km flight path with a Hexa carrying a 5000 mAh 14.8V lipo, and Astro dog tracker is about the maximum safe flight distance.
- Use the SD4000 parameters table to identify rough estimates of flight height vs. image overlap and track width requirements.
- Generally flight paths represent a tradeoff between altitude, area, and time
- It is necessary to have basic information about the changes in terrain at the site, the height of trees and other features
- Flying height based on the parameters table should be estimated from the 'highest' object in the AOI
- Current plans do not use 3D flying - but this is possible with a recent firmware update and future versions of this doc will incorporate such plans
- Consider the height of where the Hexakopter will be launched from in relationship to the tallest features
- The Hexakopter must be able to 'see' or be within line of sight of the transmitter at all times. If the hexa flies to the other side of a mountain it will likely lose transmitter signal and auto-land then crash.
Build the flight line
For the Mikrokopter, start and end the path at the launch position following a zig-zag or grid pattern, this means that little flight is wasted in merely 'travel' and all locations are imaged relatively evenly.
- Design the path so that the last waypoint overlaps / overshoots the first waypoint. Later, when the image set is trimmed to only the flight dataset there will be no gaps in photo coverage. If the last waypoint is not designed to overshoot the first waypoint, it is likely there will be a small (<10m) gap in coverage, although this does not appear to affect the point cloud.
Build the waypoint file
The waypoint file should have the format:
- FID(Default), SHAPE(Default), ID (Default), Latitude(double), Longitude(double), Radius(Int), Altitude(float), Climbrate(int), DelayTime(int), WP_Event(int), Heading(int)
- Add waypoints along the track in exactly the order the Hexa is to fly at points where the unit needs to turn
- The Hexakopter currently takes a maximum of 31 waypoints
- It is recommended to have a waypoint at least every 150m or so - to avoid excessive drift between points during flight.
- In practice, it is typically necessary to place a waypoint at each corner of the flight plan only as the Mikrokopter is severely limited in its waypoint programming. The Ardu system has no such limitation and in fact follows the path and not just the waypoints.
- Populate the ID field with incrementing integers in the order of travel - starting with 0,1,2,3,etc.
- Populate Latitude with latitude in decimal degrees
- Populate Longitude with longitude in decimal degrees
- Populate Radius with 10
- This is the Hexakopter default - we find it provides a good tradeoff between accuracy of arriving at the point and the ability for the unit to actually get to the point
- Populate Altitude with 0
- This means that the copter will not adjust altitude between waypoints.
- If a number is put in Altitude, the copter will fly to an altitude of that number plus the altitude at which altitude lock was set at.
- Populate Climbrate with 0
- This is the rate at which the copter will climb to reach the set altitude to this waypoint.
- Populate DelayTime with 1
- We find this is a good tradeoff between having the Hexakopter make sharp 90 degree turns and also having smooth travel along long transects.
- Populate WP_Event with 0
- This is used to program camera and POI functions
- Populate Heading with 0
- This is used to program camera and POI functions
- Export the table as a csv text file
- Convert the text file to a Mikrokopter wpl file using the custom script:
- python convert_ARCGIS_csv_to_mikrotool_wpl_file_v2.py waypoints.csv
Prepare the OSD image
See http://www.mikrokopter.de/ucwiki/MikroKopterTool#OSD for several options on OSD map generation
- MKTools Maps (Sometimes doesn't work and website is down)
- Open MKTools Maps
- Enter a set of decimal degree coordinates for the AOI
- Use the buttons to navigate to the area of interest and press Save Image
- Then you can open then image into the Mikrokopter Tool OSD
- GeoMapTool http://www.geomaptool.de/index.php?sid=1308497413&wx=1920&wy=1080
- Browser based.
- This seems to work pretty well, but is a bit slow.
- Just follow the on screen buttons
- OpenStreetMap can also be used to generate a background image
- Or, an image can be used from GIS
- Make a hollow rectangular AOI polygon over the image with a roughly 4:3 aspect ratio
- Make a layout with just the AOI and the image and export as JPG
- Crop the image to exactly the AOI boundary
- Load the picture in the OSD software and 'geocode' using the extents of the AOI when prompted.
- If using GCP markers - place them in the field prior to flight and record DGPS location.
- Move laptop, Hexakopter, transmitter, dog tracker, tools, cameras, batteries, spare parts to launch site.
Ground control points (GCPs)
Prior to image collection - a network of ground control points be placed within the AOI. These may be large buckets or bins. If markers will be used, collect the XYZ location of the markers with a precision GPS. If GPS equipment is not available, it is also possible to obtain the XYZ information from GIS layers, such as road or curb layers, aerial imagery, and DEMs.
Example distribution of GCPs at the Knoll site -- about 15 GCPs are distributed across the site.
- Use the orange 'Homer' buckets as GCPs
- When placing new GCPs clear/remove the grass in a 15cm dia. circle and place a steel spike flush with the ground at that location
- Spray paint the cleared area and nail with pink or orange spray paint, this will help you find it if you return in the near future.
- To save time, set up GCPs the day before by marking the locations with paint and nails and recording the GPS position
- Take a single bucket, several nails, spray paint, a GCP map and the GPS with you to the site
- At each GCP, set out the nail and paint the area, place the bucket over the nail, place the GPS on top of the bucket, and record the GPS position
- Record at least 500 positions per GCP, more is better. This equates to about 7-8 minutes of GPS time per GCP.
- Now if you come back the next day for a flight, you just set out all our buckets in the known locations and get to it!
- It is necessary to post-process the GPS positions to improve the location accuracy. Refer to the appropriate documentation for the software you are using for post-processing. Our lab currently uses Trimble Pathfinder software.
- Save this GCP list in your flight GIS so that it can be used for georeferencing and flight planning
Prepare Aircraft and Camera
- Start up laptop
- Make sure there is a microSD plugged into the Hexakopter, this is NECESSARY for recording telemetry
- Turn on dog tracker on Hexakopter and the Garmin unit
- Turn on Hexakopter with a small capacity battery that is not to be used for flight so that the unit can acquire GPS while the other equipment is being prepared.
- Check that the 3DRs are communicating
- Prepare velco straps on Hexa
- Place calibration target (Lastolite EZYBalance) in full sun
- Get camera taking photos using the simple method (a velcro strap) or the advanced method (CHDK):
- In order to use the SD4000 for image collection it is necessary to keep the camera in Continuous Shooting (CS) mode throughout the duration of the flight. Conveniently, the camera will keep whatever exposure settings are applied at the first image, maintaining constant exposure throughout the flight. We have two main methods for 'keeping the shutter button pressed' during flight:
- Simple Method
- Turn on camera in manual shooting mode
- Set camera to Program Mode (P) by switching the top switch to the middle position
- Turn the flash off (tap right on the dial)
- Adjust the function settings by clicking the center button
- Spot Light Metering is default and OK
- My Colors off
- Daylight AutoWhiteBalance is usually suitable
- ISO between 200 or ISO 400 is usually suitable
- Shutter Priority (Tv) I try to stay 'faster' than 1/800s
- Timer off
- Continuous Drive Mode
- 10MP max resolution
- Sharp Compression
- In the menu > Picture Settings:
- AF Frame: Center
- AF Frame Size: Normal
- Digital Zoom: Standard
- AF-Point Zoom: Off
- Servo AF: Off
- AF-Assist Beam: Off
- Flash off
- I-contrast: Off
- Safety Shift: Off
- Review: Off
- Disp. Overlay: off
- IS Mode: Contninuous
- Date Stamp: Off
- Set the camera to infinity focus (tap the circle left, then select the mountain)
- Lightly tape a small piece of plastic, like a small extra plastic nut from the Mikrokopter, to the shutter button (make sure not to depress the button
- Lightly wrap a velcro strap around the camera body
- Point camera at target with ISO 200 or 400 and shutter speed at least as fast as 1/800. Adjust shutter speed & ISO for good exposure on the target. Faster shutter speed is best.
- This part takes a little art: You need to use the velcro strap to push down the button while the camera is pointed at the calibration card, securing the velcro strap so that the button stays pressed. Make sure it is snug and can't wiggle off!
- Once you have the velcro strap secured, attach the camera to the Hexakopter and move on the next step.
- CHDK Method
- The 'Simple' velcro method described above is very simple, but can be challenging for some people to get right. We also suggest using CHDK (Canon Hackers Development Kit) to enabling scripting of the continuous shooting mode. This does not require using the velcro strap to keep the button pressed.
- With CHDK and the continuous shooting script installed on the SD card <maybe a tutorial on this sometime?>:
- Set SD card to lock position
- Insert SD card in camera
- Start camera in picture preview mode
- Navigate in the menu to Firmware Update and active the firmware update -- this initializes the CHDK firmware
- Load the fast intervalometer / continuous shooting script in CHDK and check that it is set to never stop taking pictures
- You need to push the Up-Left corner of the wheel to active CHDK mode -- you know this is set when additional information is displayed on the screen. This is tricky, so practice it
- Assuming other settings are set as above
- Set camera to infinity focus (left click on wheel and select the mountain)
- Point camera at target with ISO 200 or 400 and shutter speed at least as fast as 1/800. Adjust shutter speed & ISO for good exposure on the target. Faster shutter speed is best.
- Initiate the CHDK continuous shooting mode as per the onscreen instructions
- Take few photos of the Garmin GPS time for synchronization later
- Attach camera to bottom of Hexakopter using velcro straps <Need picture here of Stephen's velcro rigging method>
- FROM HERE THINGS MUST BE DONE QUICKLY AS THE FLIGHT BATTERY WILL BE PLUGGED IN AND CANNOT BE TURNED OFF WITHOUT LOSING WAYPOINTS. THE UNIT WILL REMAIN ON UNTIL THE MISSION IS OVER.
- Switch to the flight battery
- Launch Mikrokopter Tool Application
- Wait for signal to be received
- Open the OSD window
- Load the image into the OSD
- Load the waypoint editor
- Add the appropriate waypoint file
- Send waypoints to NaviCtrl
- Confirm data transfer, it should show 0/XX where XX indicates the number of programmed waypoints
- Secure dome with plastic tie
- Confirm adequate GPS status on the Mikrokopter
- At least 6 satellites are needed for flight
- When GPS status is first obtained, the kopter will beep once and the LED on the GPS will go from solid to flashing
- The OSD will read the number of satellites and a icon will appear on the image revealing the kopter location
- If you start a mission (turn the propellors on) with less than 6 satellites, no or bad telemetry will be recorded for the flight = NOT GOOD
- Time to Fly!
Fly the Mission
- Move Hexakopter to open launch area
- Make sure unit is level by placing bubble level on one of the arms
- Point red arm in direction of flight lines use OSD compass icon for reference
- Set 3-way switch to OFF (all the way down)
- Set Altitude lock switch to OFF (down)
- Calibrate units level by moving left stick up and to the left until beeps
- Return left stick to bottom center
- Turn on props by moving left stick to lower left
- Take off
- At 20m or so in altitude turn on GPS lock to keep horizontal position by moving 3-way switch to the middle position
- Climb to desired altitude - read off of both OSD altitude and a laser range finder for backup - rely on laser for the most accurate reading
- Set Altitude lock by moving Altitude switch to the up position
- When ready set unit to follow waypoints by moving 3-way switch to top position (marked with HOME sticker and closest to flyer
- Watch the unit fly its path
- When unit has finished its course and has stopped at its last position tell it to come home by quickly moving 3-way switch all the way down (OFF) then all the way up (HOME)
- The unit should now move to a hover above the starting point
- Set the 3-way switch to the middle position to keep horizontal control
- Slowly move the left stick down until unit starts to descend
- At a comfortable altitude (about 30m) move left stick to the middle and turn off altitude lock
- If altitude lock is turned off while the left stick is down it well descend too rapidly
- Slowly manuever the unit down by gradually moving the left stick down
- At about 20m or so release the position hold by moving the 3-way switch all the way down (OFF)
- Bring the unit to a soft landing
- Turn off Hexakopter
- Turn off transmitter
- Turn off camera
- If this is the last flight
- Turn off dog tracker & Garmin
- Transfer photos to the laptop
- Remove dome and remove microSD
- Transfer GPS telemetry to a job folder
- Replace MicroSD card
- Turn off laptop
- Pack equipment, ground control markers and head home
Extract Images and Telemetry
- Use telemconvert.py script to convert the telemtry GPX to a csv file
- python telemconvert.py GPS.gpx
- Load CSV file into ArcGIS and display XY in WGS84 coordinate system
- Note - I have found that if you export the CSV to a shapefile now the time field is corrupted and lost by ArcGIS - the only solution I have found is to try to use a geodatabase, otherwise keep the imported CSV XY points as an 'events' layer
- Using the time stamp on points and the camera time read with the EXIF Reader software identify the start and end points of the flight
- First it will be necessary to calibrate the camera and GPS time -refer to pictures of the GPS time for this difference.
- It may be helpful to symbolize the track by NCFLAG to note where the unit arrived at a waypoint: Code 24 means 'arrived at waypoint'
- One useful technique may be to select a subset of photos that includes the first 'turn' of the copter along its path from the start and then process this small set quickly in the computer vision software. When the model is complete, you can then see which photo was located at the 'turn' and can then synch the timestamp of that photo with the time stamp of the GPS point at that 'turn'.
- Segment out the flight photos from the takeoff and landing photos based on time
- Process photos within Photoscan to generate a point cloud
- Refer to the georeferencing document for instructions on georeferencing and noise removal of the point cloud data.
Generate 3D Point Clouds with Computer Vision
Run Computer Vision Software
Once the photos have been trimmed to the flight set it is time to make point clouds!
These instructions are for running Agisoft Photoscan v0.8.5 build 1423 (25 April 2012). We often use a screen capture software like Cam Studio to monitor the status of the job by placing the progress indicator in a small area near the system date/clock and making a screen capture every 30 seconds. This makes it possible to identify when the job actually finished.
- Open Photoscan
- Go to Workflow > Add Photos
- Add the flight photos
- Go to Workflow > Align Photos
- Use Accuracy: High, Pair Selection: Generic
- Click OK and let it run!
- Once Photoscan has completed, you will need to SAVE the job or else lose the work that was done
Running the structure from motion software Bundler requires that you first compile the program on your system. This tutorial does not yet cover the details of compiling Bundler on your system. However, once Bundler is compiled the process is fairly simple.
- Execute the RunBundler.sh script from within the directory of your images
- We made a slight modification to the script to produce a time stamp in a text at each major stage of processing, to provide a way to keep track of how long the computation takes.
- Once Bundler is completed it automatically generates the PLY 3D point cloud file that will be used in further computation.
Export 3D Point Cloud and Camera Points
Exporting Data from Photoscan
To export the point cloud:
- Go to File > Export Points: and select a PLY file for export.
- Generate an ASCII PLY file, not a binary file, and export colors, but not normals
To export the cameras:
- Go to Tools > Export Cameras
- Generate an XML file or other file type, for example CHAN.
- In the next steps we will generate a camera list from this output for further work
Export Data from Bundler
By default, Bundler will produce the necessary output upon completing the script, however, you can generate a list of just the camera points easily.
- Load the point cloud in GIS, or Excel, or as a Numpy array in Python
- Camera points in Bundler output will satisfy these condition based on color:
- Red: red=255 green=0 blue=0
- Green: red=0 green=255 blue =0
- Yellow: red=255 green=255 blue=0
- The actual camera centers are the Red and Green points. For every red or green point, there is a yellow point that defines the camera's orientation.
Georeferencing Point Clouds
Before the point cloud can be used for analysis - it is necessary to perform a georeferencing and noise removal on the dataset.
This process involves two main steps:
- Generating a table of 'coordinate pairs' to base the georeferencing on.
- Executing a custom python script to sequentially georeference the whole point cloud and perform noise removal
A note: Here, ground truth coordinates are described in WGS84 UTM meters and code has only been tested for this coordinate system.
Ground control points (GCPs)
- It is recommended that prior to image collection - a network of ground control points be placed within the AOI as possible. These may be large buckets or bins. If markers will be used, collect the XYZ location of the markers with a precision GPS.
- If GPS equipment is not available, it is also possible to obtain the XYZ information from GIS layers, such as road or curb layers, aerial imagery, and DEMs.
Prepare for Georeferencing
- Measure Buckets
- Extract GPS Lists
- Use the telemconvert.py
- Need to convert GPS coordinates to UTM
- Convert MK pressure-based altitude readings to Elevation
- Based on our own measurements: Y (meters) = 0.0393(Altimeter) + 1.3214 + (Elevation of flight start location)
- Extract CAM Lists
- From Bundler
- Segment last points, take the red and green points only
- From Photoscan
- Export as XML and run our parse_pscan_xml_for_cam_location_v2.py code
- Export as CHAN and take the first 3 columns
- From Bundler
Ground Control GCP 'Bucket' Method
- Have created a 3D point cloud using computer vision software and exported as PLY
- If using Photoscan - convert the PLY into ASCII format using Meshlab by opening the PLY in Meshlab and choosing Export > PLY > and uncheck binary encoding
- Have installed Menci Scanview software http://www.menci.com/zscan/index.php?option=com_content&task=view&id=45&Itemid=26
- Open the ASCII PLY in a text editor and remove the header rows. Save as a new TXT file
- In Menci Scanview, Import the new TXT file
- Next, prepare a GIS point shapefile or an Excel table for adding in the coordinate pair GCPs. It should have a format like this:
Point# CloudX CloudY CloudZ UTMX UTMY UTMZ 1 -11.064 -28.094 20.266 435294.1 4428618 66 2 14.809 -3.086 28.378 435359 4428638 66 3 -6.427 36.571 30.679 435409.2 4428569 62 4 -39.569 2.011 22.785 435318.7 4428544 65
- Where CloudX,Y, and Z are fields for the coordinates of the ground control points in the point cloud coordinate system and UTMX, Y, and Z as the ground truth coordinates in meters of the ground control point - as measured with GIS or DGPS.
Measuring point coordinates in Scanview
Scanview is a 3D model visualization tool that also contains a model measuring tool. This makes it easy to obtain the point cloud coordinates of GCPs based on 3D structure and color. But first, a quick note on navigating in Scanview:
Command Function Left Click and hold Rotate the model in 3D Right click and hold Zoom Middle Click and hold Pan Ctrl + Left click and hold Rotate current display in 2D (like turning the whole monitor) Ctrl + Shift + Left Click Center the origin of rotation at the click point
- Start by navigating the point cloud to a recognizable viewpoint (this will almost certainly require all of the above commands
- Open the Measure tool (the little ruler in the toolbar)
- Navigate to the location of a GCP (either a distinct feature or a placed marker)
- Shift + Left Click to place a measurement at that location
- It is possible to make many measurements and average them, or just use one.
- Record the point cloud XYZ coordinates of the GCP in the table or shapefile along with the ground truth UTM coordinates.
- Repeat for the rest of the ground control points
- It is recommended that you generate a list of at least 12-15 GCPs, distributed across the study area. Only 5-6 GCPs will be used in the georeferencing and the rest will be used to evaluate the georeferencing results.
When this is complete, you should have a list of the point cloud coordinates for the GCPs and corresponding coordinates from the GPS. Make this list into a simple, tab-delimited text file:
- -11.064 -28.094 20.266 435294.1 4428618.0 66.0
- 14.809 -3.086 28.378 435359.0 4428638.0 67.0
- ... ... ... ... ... ...
Then execute the AutoHelmert script:
- python AutoHelmert_v1_5.py GCP_list.txt <optional 7-parameter initialization file>
By default the script uses a generic set of parameters to initialize the optimization (rotations: 1, scale: 1, translations: 350000, 4340000, 100). If you would like to provide an alternate parameter file, for example to attempt to improve the optimization, you can do so by adding in the name of the text file at the command prompt.
This script will then perform 4 steps:
- Find the center most point and five additional points that represent the furthest points from the center
- This is done by finding the combination of points that have the greatest cumulative length from the center point
- Perform a least-squares optimization of the Helmert 7-parameter coordinate transformations to convert point cloud coordinates to geographic coordinates
- Evaluate the error of the solution based on the remaining GCP markers that were not used in step 1
- Export the list of optimal parameters, and several text file and graphical reports of the solution, including RMSE and NSSDA error statistics
The parameter file can then be used in the Point Cloud Filtering step below.
Aerial Path Spline Method
The spline method for georeferencing was developed to enable geocorrection of Ecosynth point clouds when it is not possible to place GCP markers at the site or when GCPs would not be visible, (e.g., in a closed canopy forest).
The idea behind this technique is that the GPS path from the Hexakopter telemetry (trimmed to the actual flight) and the path represented by the sequence of camera/photo points from the computer vision software should have the same relative geometry and that the 7-parameter transformation can be solved based by fitting the two paths together.
To accomplish this, the entire aerial GPS track (UTM coordinates) and the entire set of camera positions along the flight path (SfM coordinate system) are fitted to independent spline curves, from which a series of 100 XYZ pseudo-pairs of GPS and SfM camera locations are obtained using an interpolation algorithm (Python v2.7.2, Scipy v0.10.1 Interpolate module) and then used as input for the georeferencing of point clouds using the same Helmert transformation algorithm used in the GCP method.
To perform georeferencing with the spline method:
- Reduce the GPS path to a simple, tab-delimited text file of format, where the column values are UTMX, UTMY and elevation Z in meters, do not use a header:
- 435294.1 4428618.0 66.0
- 435359.0 4428638.0 67.0
- ... ... ...
- Place this new GPS_file.txt into a new folder
- Place the XYZ CAM_list.txt into the same folder
- Run the spline command:
- python splines.py GPS_list.txt CAM_list.txt
The result will be a report of the internal horizontal and vertical RMSE error of the least-squares optimization as well as a text file of the solved 7 parameters. This 7-parameters info is the most important output of the georeferencing step and is used to georeference the entire point cloud in the next stage, Point Cloud Filtering.
Point Cloud Filtering
We are almost done! Before the data can be processed into a CHM, DTM or any other product, it is necessary to apply a final stage of processing and filtering to the point cloud. This stage performs the following tasks and generates the following output:
Function Parameters Output Coordinate transform of point cloud 7-parameter file
raw ASCII PLY point cloud file
Exports unfiltered, but geocorrected point cloud Sort out camera points looks for Bundler cameras by RGB value Exports a set of camera points and a set of non-camera points Global X & Y Filter
(Crop to AOI)
Zscore/SD cutoff value (default:3)
(AOI XY Extents)
Internal point cloud file where extreme X and Y values have been trimmed away
(A soon to be updated version of the script uses an AOI (Area of Interest) as input for cropping the point cloud)
Local Elevation Zscore filter Zscore/SD cutoff value (default:3)
local window size (default:10m)
Exports a new point cloud where height values within 10m x 10m grid with |Zscore| > 3 are discarded Local Median Filter FilterType (default: median)
local window size (default:1m)
Exports a set of points matching the window size, where only the median elevation value is retained, XY location is the window center
This process is executed entirely within the Python noise_removal.py script. The script primarily utilizes functions from Numpy and Scipy to perform manipulations, queries, and statistical computation on the XYZ-RGB point cloud as an array dataset. Currently (2012-09-08) it has some memory issues with large point clouds when run in 32-bit in Windows. This is because it was programmed a little lazily and runs just fine in 64-bit implementations of Python in Linux.
To execute the script:
- python noise_removal.py point_cloud.ply params.txt
All output is numbered and named based on the stage of the program and located in an 'output' folder within the current working directory. When the script is run successfully, the result is two point clouds (local filter and median filter) that have been georeferenced, noise filtered and are ready for GIS (with the addition of a header, thanks ESRI!) or further analysis. In addition, the script generates PLY files of the output point clouds so that the results can be immediately previewed in Meshlab or a similar software.
Generate a Digital Terrain Model
A digital terrain model (DTM) is a raster data layer that represents the ground surface elevation, typically when aboveground features (e.g., trees, forest canopy, and buildings) have been removed through filtering of 3D point clouds. This filtering can be performed with many different tools and algorithms, but often requires as much art as science in order to get the tool/algorithm parameters to work effectively for a dataset.
A DTM is a necessary component of a 3D remote sensing pipeline for measuring tree height. The elevation of the terrain is subtracted from the elevation of canopy points in order to determine true height above ground.
Previously, we used the ALDPAT LIDAR point cloud filtering toolset to apply terrain filtering algorithms to LIDAR and Ecosynth point cloud datasets. Recently, other tools have become available, e.g., LASTools and may be more effective and automated for terrain filtering.
To use a terrain filtering algorithm, import the median filter point cloud text file from the point cloud processing/noise removal output folder and execute the algorithm. A good result will maintain a large number of points and a high mean density (e.g., no large gaps in coverage) and have minimum error compared to a high quality standard like a LIDAR DTM.
Once a point cloud has been filtered to include only ground points, it is necessary to generat a raster DTM.
- Add a header to the output file if there is not one already -- typically X,Y,Z is sufficient
- Import the layer into ArcGIS, create XY points, and generate a new shapefile of the ground/terrain points
- Use interpolation, e.g., Kriging, to generate a raster model from the ground points layer
- Crop the raster model to the project AOI
In addition to generating a DTM using a filtering algorithm, it is also possible to use an existing LIDAR DTM. When provided by a LIDAR contractor, these are typically called bare earth products because the contractor has already filtering off the above-ground points.
Generate a Canopy Height Model
A canopy height model (CHM) can be a raster or point cloud (vector) dataset where elevation Z values represent the height of that pixel/point above ground. This is different from the output of the point cloud filtering process that generates a point cloud representing the surface model where each elevation Z value represents height above mean sea level (or geoid or ellipsoid, depending upon the reference datum used).
The CHM is an important data product of the Ecosynth 3D remote sensing pipeline as it is a common dataset used for estimating biomass density, parameterizing fire models, mapping tree crowns and gaps, and predicting species diversity.
Creating a CHM is very simple once a DTM has been created or if an existing LIDAR DTM/bare earth data product is used.
- Add a header to the the local filtered point cloud text file that is in the output folder of the point cloud filtering process
- The header can be simple: X,Y,Z,R,G,B -- but is necessary for ArcGIS to be able to import the file
- Import the file into ArcGIS as a text file, and create XY events using the X, Y, and Z fields and in the correct coordinate system. Refer to the ArcGIS help if you are unsure about importing XY coordinates.
- Export the XY events layer as a new shapefile
- Add in your DTM raster model for this area
- Use the Spatial Analyst tool 'Extract values to points' to assign a DTM elevation value to each point cloud point -- this creates a new shapefile with a field RASTERVALU that represents the DTM pixel value below that point
- In the new layer that is created from the last step, add a new field 'HEIGHT' (float) and subtract the RASTERVALU from the point Z value
- Ta-da! You have now created a layer that shows canopy height instead of elevation!
From here you can perform several additional tasks, including creating a CHM or extracting point cloud height statistics per forestry plots, as we did in our papers.
Compute Point Cloud Grid Stats
If you want to generate a raster CHM from this dataset, it is recommended that you perform an additional stage of analysis to remove some of the complexity of the data. In this stage the CHM point cloud will be analyzed to extract information about point cloud statistics (count, density, mean/max/min/etc height) within 1m x 1m bins. This is useful because kriging or other types of analysis of the millions of points in a point cloud is very difficult for ArcGIS. By reducing the number of points in the point cloud, it is a lot easier for ArcGIS to produce useful raster maps: e.g., point density maps, CHM maps, etc.)
To extract grid statistics from the point cloud:
- Export the CHM point cloud table as a text file
- From the command line, execute the grid_stats script:
- python grid_stats_v3.py CHM_point_cloud.txt
The result will be a text file point cloud with X and Y values representing the center of 1m x 1m bins across the extents of the study area, within point cloud statistics for each bin. This can then be added back into ArcGIS as an XY layer, which can then be exported as a shapefile and interpolated into CHM and density maps.