Lidar analysis

Introduction

Summit elevations have traditionally been taken from one of two sources.  The first is topographic maps, where elevations were often derived from stereo photography, or much more rarely, from an historical on-site survey.  The second, more recent source is digital elevation maps (DEMs), which are regular grids of elevations.  These are typically collected from satellites or aircraft and based on radar or other remote sensing.  These methods have several significant sources of error:

• The inherent error in topographic maps is often large, on the order of one contour, which can be in the tens of meters.

• Topographic maps are often many decades old, meaning they don't capture changes to summits to due mining, volcanism, etc.  The sea level reference used on older maps can differ by several meters from the most recent reference, depending on location.

• DEMs are usually surface models, meaning they include treetops, buildings, etc.

• DEMs with broad coverage usually have large spacing, such as 30 meters, meaning they tend to underestimate summit heights, sometimes significantly.

Data collection by Lidar attempts to address all of these problems.  It consists of recent measurements with very high density and inherent precision.  At infrared wavelengths and high point density, it can often penetrate vegetation to reach the bare ground.  Because of its cost, it is only available in a small number of countries.  But where it exists, it can dramatically improve summit elevation estimates.

Computing prominence from Lidar

In 2022 I extended my existing prominence code to take the USGS 1m DEM as an input.  The main difficulty here was that the 1m DEM was in a different projection, namely Universal Transverse Mercator (UTM).  Previously the code had worked only on data in the lat-lng projection, which was in very common use for datasets with large coverage.  But UTM is also a widely-used projection, so I figured it would be worth the effort to support it directly in the code.   I also sped the code up where I could, as the 1m DEM is 100 times larger then the 10m DEM I'd previously worked with in the US.  I analyzed a few small areas in Arizona, which made it clear that there were new undiscovered P300s (peaks with 300 feet of prominence) in the data.  However, 1m data wasn't too widely available, so there things sat for some time.

In early 2024 I thought it was time to take Lidar analysis more seriously.  Data was becoming more widely available, obsoleting analysis done by simply looking at topographic maps.  There was one glaring obstacle to doing this at large scale.  Even within a single delivery of data from the USGS, different areas were in different projections (called different UTM "zones" in the lingo).  In Europe, each country was using its own specialized projection.  The were clearly too many projections for me to build support for all of them into my code.  Instead, I wrote scripts that would take data in an arbitrary projection and convert it to a single common projection that my code supported.  The lat-lng projection was the obvious target.  The reprojection itself is done by an open-source program called GDAL.  This process introduces some error, because the original height samples are combined with their neighbors during the reprojection.  As will be seen, this was a relatively minor source of error.

If you had Lidar data for the whole world and you were building a new peak database from nothing, that would be basically the whole story: just process all the data, and you'd get all the mountains in the world with high accuracy.  But in reality, Lidar data is collected by a wide array of governments and other agencies, resulting in boundaries that are often inconvenient, and sometimes hard to fathom.  For example, here are the boundaries of Lidar datasets for the majority of Utah:

In order to compute the prominence of a peak, the elevations of both the summit and the key saddle must be known.  If the summit is within a dataset but the key saddle is not, the computed prominence will be wrong (in particular, it will be too large).  Of course, it's not possible to know a priori whether the key saddle is within a given dataset.  One can use an existing peak database as a hint of the saddle location, but not as a guarantee.  After all, one of the purposes of the Lidar analysis is to correct errors in key saddle locations.  A further complication is that I was applying my analysis results to the existing Peakbagger database, which contained over 100,000 peaks.   If a peak had two summits A and B of nearly the same height, where A was in Peakbagger, but my analysis showed that B was higher, it was not enough to simply add B to the database.  Peak A must also be adjusted.  This was not always possible to automate, since A might have very tiny prominence, less than the minimum cutoff used in my analysis.

One implication is that results are better when the area being analyzed is large, so that more peaks will have their key saddles within the dataset.  In practice, this means waiting until large contiguous areas are available.  The other implication is that the prominence results require a tedious manual review.  Inside a single window, I loaded a topo map, all existing Peakbagger peaks, and all P300s identified by the prominence analysis.  I paid special attention to several particular cases:

• Existing Peakbagger peaks where the elevation or prominence changed dramatically.  This could be a due to several reasons, from a peak's elevation being entered incorrectly, to a key saddle that was outside the dataset, to a peak that had been lowered by mining.

• Existing Peakbagger peaks where the key saddle moved more than a minimum distance.  These were usually inaccuracies in Peakbagger, but sometimes they were a sign that the key saddle was outside the dataset.

• New P300s that were close to P300s in my previous global analysis.  These were usually valid, but I still checked each one by hand.

• New P300s that were not close to a P300 in my previous global analysis.  These were often invalid, usually because their key saddle was outside the dataset, but sometimes they were valid, and simply differed in location from the global analysis due to the Lidar data being more accurate.

For any new peaks that I found, I wanted to use the best possible name, rather than defaulting to simply using the elevation as the name.  In the U.S., I cross-referenced the each new peak with official names from the Geographic Names Information System, and the names of benchmarks from the National Geodetic Survey.  Even so, I still checked the topographic maps for additional names.  There were names on the map not present in the official databases, and names in the databases not on the map.

While the USGS has plans to publish a 1m DEM for the entire U.S. (outside Alaska), today there are still many gaps in coverage.  In a few cases I was able to fill holes with Lidar data available from states or even counties, but the cost in additional time was high.  These smaller agencies had wildly differing projections, units (feet versus meters) and ways of accessing their data that varied from the baroque to the broken.  It would sometimes take months of back-and-forth over email to decode how to download their data.  However, it was worth it to cover critical missing areas of Colorado, New York, and New Hampshire, where otherwise peaks and prominence values would have gone missing.

Countries in Europe with Lidar available presented their own challenges.  Just getting the data on my local machine was sometimes a struggle.  Switzerland was covered by 0.5m data in over 43,000 tiny files totaling over 700 GB, while Slovakia is covered by a single file over 150 GB in size.  In Germany, each state collects and hosts its own data, using completely different systems, naturally.  And France---well, let's not talk about France.  Strangest of all was Austria, where the datum (the sea level reference surface) had changed since its official topographic maps were last updated, and thus the Lidar elevations had to be adjusted according to a special recipe.

Finding the ground

In an ideal world, this would be enough.  But I soon found that some peak elevations were coming up suspiciously short of known elevations.  I understood that the DEM was generated by averaging or otherwise combining the raw points measured by the laser within a 1m square, which could result in a small loss of elevation at a summit.  But the errors I was seeing were far larger than this.  It was time to go deeper and look at the raw "point cloud" generated by the Lidar device.

To understand what I found, it's necessary to describe more details about how Lidar data is collected and processed.  (This is still a simplified picture.)  A plane fires an infrared laser at high frequency downward, and records the strength and time delay of the return.   From this it is possible to calculate an intensity value and an elevation for a point below.  The laser may have bounced off a building, a rock, bare earth, or a tree.  Usually multiple passes are made, making it possible to "see" different sides of a building, or both the leaves of a tree and the ground below.  The set of points might look like this:

The accuracy of the calculated elevations is generally excellent, on the order of centimeters.  But now comes the crucial and most complicated step.  Each point must be "classified" (labeled) as either ground, vegetation, water, building, or perhaps "unknown" if the classification algorithm isn't able to figure out what the point is.  Classification of points as ground or non-ground is mostly independent of the Lidar hardware, though of course results will be generally be better the denser the cloud of points.  It turns out that doing this well is very difficult, and in fact it is still an unsolved research problem in general.

When the point cloud is converted to a regular one-meter grid (the DEM), only ground points are included.  Therefore, if some ground points are incorrectly classified as non-ground, the elevation in the DEM can be too low.  If some non-ground points are classified as ground---for example, the leaves of a tree---the DEM will be too high.

As it turns out, the algorithms used by the companies collecting Lidar data for the USGS are very conservative; that is, they try very hard not to accidentally interpret a tree or a building as ground.  In fact, they are so conservative that they occasionally fail to classify the very top of a mountain at all, meaning that the DEM at summits is too low.  This is the cause of the low elevations I was seeing.

The problem turns out to be both wide and deep.  An analysis by the ListsOfJohn team showed that the median underestimation above treeline on a set of peaks in Colorado was over 1 meter---hardly the centimeter-scale accuracy one  hopes for.  They also found some cases where classification failed on an entire ridge, underestimating the summit by dozens of feet.

Left: Lidar point cloud for a treeless ridge in Wyoming.  Center: Ground points according to the USGS classifier.  Right: Ground points according to a different classifier.

With further research, it might be possible to come up with a better classification algorithm and build new DEMs that more accurately extract summit elevations from the point cloud.  Alternatively, one could manually review every summit by hand and attempt to identify the highest ground point by eye.  Both of these approaches are prohibitive due to the enormous size of point cloud files, and the difficulty that even experienced humans have in telling rocks from trees in Lidar data.  Clearly, some kind of compromise is needed.

I found a point cloud classifier with a different algorithm, one that is far less conservative.  In fact, it is not conservative enough: it often classifies trees and buildings as ground, leading to summit elevations that are too high.  Thus, it is only useful above treeline or on very rocky peaks without vegatation.  There must also be no buildings or antennas on the summit.  Determining whether this classifier is safe to use thus relies on manual review of satellite imagery.  It might be worth attempting to automate that part in the future.

I built a new semi-automated pipeline that works like this:

1. For each peak, use an index file provided by the USGS to determine the point cloud file containing the peak.

2. Download the point cloud file.

3. For better performance, crop the point cloud file to a small radius around the peak (say, 50 meters).

4. Run the new classifier on the cropped region.

5. Take the highest classified ground point and output that elevation and position, along with an identifier for the peak.

6. Manually review satellite imagery of the summit to determine whether it is safe to accept the new elevation.

Strangely enough, the most difficult step is the first.  The index files provided by the USGS appear to have no standard schema.  They are often missing, misplaced, misnamed, or broken in other creative ways.  There are still small areas where I haven't been able to run this pipeline due to a lack of a working index file.

While this manual review has been able to identify and compensate for some pretty egregious errors in the DEM, it is unsatisfying and tedious.  Another weakness is that it adjusts summit elevations after prominence has been computed.  It is possible that adjusting a summit elevation in this way could cause a peak to become taller or shorter than a neighbor, causing them to swap key saddles, but I don't have a way to detect this.  A further deficiency is that I don't review manually saddles for errors.

Future work

The most serious drawback of existing Lidar data is incomplete or erroneous ground classification.  This is an active area of research.

Lidar point densities are increasing with time.  Older data can have densities as low as 1 point per square meter, which can miss summit boulders in some important cases of interest, such as Challenger Peak, a potential high point of Alameda County, California.  Newer data has much higher densities, 10 points / square meter or more.  It may be worth reanalyzing certain areas in the future.

Published Lidar coverage in the US (left) and Europe (right) as of September, 2024