Introduction

We describe here how the Jura Mountains mapping website uses GEDI LIDAR (light detection and ranging) transects (PDF).

But first, some guidance for users.

Images giving the profiles of the density of biomass  - technically the Plant Area Volume Density (PAVD)  - along transects are made available by selecting "Transects" in the map menu. The map should then be rotated (use Alt + Shift keys) and zoomed to the resolution that makes the length of the transect (the blue line below the transects PAVD image) equal to the PAVD profile shown in the image.


Project level

OpenStreetMap technology is rarely applied at the local project-level as opposed to the global level.  The Jura Mountains mapping website aims to overcome this shortcoming by piloting open-source OpenStreetMap project-level applications.

A recent development for application at the global level involves the  Global Ecosystem Dynamics Investigation (GEDI) high-resolution earth-observation LIDAR data.

GEDI Level 1 and Level 2 products have been made available since late-2019. They allow the canopy height and vegetation density of forested areas to be estimated.

Level 3 products that grid spatially interpolated Level 2 footprint estimates of canopy cover, canopy height, leaf-area, and vertical foliage profiles were planned for release in mid-2020, and are presumably still planned for release.

Level 4 products - the highest GEDI product level - will eventually use calibrated models to derive estimates of above-ground biomass density footprints.

These footprints will then be used to estimate biomass within 1 km square cells in order to monitor the Earth's carbon cycle at the global level as a component of carbon accounting aimed at reducing the level of atmospheric carbon dioxide which drives climate change.

The three GEDI lasers mounted on the International Space Station fire 242 times each second to produce eight beam tracks which are spaced approximately 600 m apart in the cross-track direction (see specifications). Two beams are full-power beams and a third laser beam is split into two "coverage" beams. Every other shot of each of the beams is optically dithered across-track to give an additional beam tracks with irregularly spaced shots. Each beam samples the Earth at points that are approximately 60m apart along the beam track's orbit, resulting in almost continuous along-track coverage. Each illuminated footprint reflects LIDAR energy with a diameter of about 30 m and a vertical measurement accuracy of about 2o to 3o mm.

The focus at the project-level pilot is to evaluate the application of GEDI products to the built environment that comprises urban and nature-based infrastructure which blend vegetation and man-made structures.

As a first step, Level 2 GEDI metrics for transects through the Jura Mountain forests are being examined. The aim is to see if the vertical foliage profile can be used not only to monitor logging and other forestry activity but also to spatially quantify the built environment. 

Web mapping uses OpenLayers since map rotation is needed to align images of PAVD profiles with GEDI transects. Drop-down menus are used now to select profiles, but this will be changed at some stage by using pop-up images.


Implementation

Implementation is fairly straightforward. As regards data processing, the starting point is to map vegetation profiles along sections of the GEDI beam tracks that span forests in the Jura Mountains.

Images of PAVD profiles are obtained using a Jupyter Notebook published by the US Land Processes Distributed Active Archive Center (LP DACC). 

As the Notebook explains, using a terminal (search for "cmd" , the Command Prompt in Windows 10), create a Conda environment and install and run the Notebook:

  • conda create -n geditutorial -c conda-forge --yes python=3.8 h5py shapely geopandas pandas geoviews holoviews
  • conda activate geditutorial
  • conda install jupyter notebook
  • jupyter notebook

Installing Conda can be a bit troublesome. For Windows 10 with Python already installed, Miniconda is adequate. After downloading and following instructions, let the installer set the PATH variable, which can checked with the terminal command:

  • echo %PATH%

GEDI Level 2B data (as HDF5 formated files) then needs to be downloaded. The DAAC tutorial Notebook demonstration file is:  GEDI02_B_2019170155833_O02932_T02267_02_001_01.h5

As explained in the Notebook, a NASA Earthdata  account is needed and the .h5 file must be downloaded or copied into the same directory as the Notebook (for a simple installation, the Notebook and the .h5 file will be in c:/Users/admin where admin is a Windows 10 administrator account).


The Notebook

In running the Notebook, an error will probably be thrown when the PAVD profile of a single shot is printed, so the single-shot PAVD profile is not displayed - but it can be saved as a .png image by adding to the relevant Notebook cell:

  • hv.sav(path1,'path1.png') 

Second, care must be taken in specifying the GeoJSON file that is used to define the boundary of the beam track that will be sampled. The tutorial explains that a GeoJSON file for the Redwood National Park must be downloaded from Earthdata (or other sources - see below) to the same directory as the Notebook and the .h5 file. 

The Notebook in fact uses the maximum and minimum dimensions of an envelope that encapsulates the GeoJSON. It is usually wise to define the sampled area as a simple polygon. However, the GeoJSON file must be a denoted as a MultiPolygon. Fortunately the GeoJSON can be expressed in Web Mercator EPSG 4326 units (i.e., in terms of coordinates given by say Google maps).

Once the GeoJSON is loaded in the Notebook, the full GEDI orbit is plotted. By zooming the orbit to the area outlined by the GeoJSON file one can select the beam shot that will become the centre of the PAVD profile. Clearly, selecting a point that is at the centre of the intersection between the GeoJSON polygon and the orbit is the most sensible approach. However, by default (which can be changed), the PAVD profile is generated for a scan that covers 50 laser shots on either side of the selected shot. So different parts of the demarcated area can be sample by selecting different shots. The selected beam shot must then be entered in the Notebook.

The last section of the Notebook exports the section of the GEDI beam swathe that was used to calculate the PAVD profile across the area designated by the GeoJSON and centred on the selected shot. The swathe has (as shown on the Jura mountains GEDI map) two full-power beams with a third laser beam split into two coverage beams. Each of the beams is optically dithered across-track to give four additional beam tracks. The beam number for one of the eight beams - usually a full-power beam - must be entered in the Notebook.


Some observations

The Plant Area Volume Density (PAVD) - the derivative of the cumulative Plant Area Index - that measures the volume density distribution of biomass is recognised as a sensitive measure of forest dynamics.

In the figures, the PAVD is plotted against height in order to show the vertical distribution of biomass in forested areas. Several strata with high PAVD values are generally of interest, namely biomass concentrated as follows:

  • bottom - understory (trunks, shrubs and juvenile trees); 
  • middle - sub-canopy (branches and leaves of sub-canopy trees and trunks of the canopy trees in the understory); 
  • top - tree canopy (canopy branches and leaves).
While GEDI's LIDAR can penetrate forest canopies, optical imagery certainly cannot. So confirming conclusions drawn from PAVD observations generally requires ground-truth. With this in mind it is noted that for the transects displayed so far, which are for predominantly coniferous forests in leaf-off mid-winter with some snow cover at higher altitudes:

  • high top, middle or bottom densities can extend over a considerable distance (200-300 m);
  • immature forests have densities that vary litle with height;
  • some mountain tops appear to have exceptionally dense sub-canopies (snow-laden trees?).

Ground truth will establish whether or not these features, especially the first, reflect forestry practice.

The impact of topography on PAVD is well known, so it is probably not surprising that canopy heights are underestimated considerably for very steep slopes.  One assumes that this means that the vertical distribution of biomass appears concentrated into a small height. At what point this effect may impact the identification of man-made structures is unclear.

This concern reflects a more general issue relating to high-resolution LIDAR earth observation of mixtures of man-made structures and vegetation that is discussed below. 

Jura Mountain forested and non-forested areas (generally mountain pastures with scattered trees) are not unsurprisingly easily resolved. GEDI transects also clearly resolve man-made features down to the scale of 2m wide earthen tracks that are not tree covered. 


GEDI data

The best way to identify GEDI orbits that intersect the area of interest (technically in terms of the Jupyter Notebook, the input bounding box) is to use the GEDI Finder. In the case of the Jura Mountains:

  • https://lpdaacsvc.cr.usgs.gov/services/gedifinder?product=GEDI02_B&version=001&bbox=[46.35,6.0,46.75,6.5]

The response is a list of linked files equivalent to the link for the .h5 file used for the Notebook:

  • https://e4ftl01.cr.usgs.gov/GEDI/GEDI01_B.001/Year.Month.Day/GEDI01_B_2019219181545_O03694_T02694_02_003_01.h5
GEDI L1, 2A and 2B products are also available as follows:

  1. ESAMAAP - under the GROUND DATA dropdown;
  2. LP DAAC - as noted in the GEDI Quick Guide (PDF)
  3. Earthdata

Publically available GEDI data started to be available in April 2019 and the latest data is for January 2021.

GEDI Pull is a Python wrapper for GEDI Finder that automates a combined searching and downloading  of GEDI data from Earthdata.

The GEDI Algorithm Theoretical Basis Documents (ATBD) provide the basis for GEDI products. Level 4 documents and products covering gridded metrics have yet to be released. Product leads are listed.

DAAC has released Version 1 of  Level 3 products (Gridded Land Surface Metrics) but the link to the ATBD is missing although the ATBD is available in the DACC data package (logon required). The products give the mean canopy height, standard deviation of canopy height, mean ground elevation, standard deviation of ground elevation and counts of laser footprints per 1 km x 1 km grid cells. This first release of L3 products are simple averages and standard deviations of the footprint profile metrics within each 1 sq. km cell. Data with a user-selected projection and format can downloaded using a data access tool.


Developments

GEDI has the highest resolution and densest sampling of any LIDAR instrument in orbit to date. It has a nominal mission length of two years, starting in early-2019. The GEDI Level 2A data product demonstrates the enormous amount of data that is generated. The product is based on 156 layers for each of the eight beams, and includes ground elevation, canopy top height, relative return energy metrics (describing canopy vertical structure, for example), and many other interpreted products from the return waveforms.

Applications to date have tended to aim to improve maps of forest height and biomass across the globe. A notable derived service is Global Forest Canopy Height that can be downloaded. A programme to establish calibration plots has been set up. It seeks coincident small-footprint (at least 25 m in diameter) LIDAR data and field-measured plot inventory data. GEDI's success in measuring globally the mean biomass in 1 km cells will probably depend on the outcome of this calibration programme. Technologies based on Entwine and Caesium data storage and display and LOPoCS point-cloud streaming from pgpointcloud-extended PostgreSQL databases to serve biomass mapping will also be crucial.

With regard to the application of OpenStreetMap technology at the local project level, the main issue is whether or not GEDI and its successors will help in spatially quantifying the built environment in terms of vegetation and surface characteristics. The hypothesis here is that the construction and operation of an infrastructure project should be assessed, at least in part, in terms of the project's impact on its microclimate. 

Several approaches for classifying vegetation and surface characteristics are used, notably GIT, LULC, LCZ, HERCULES and USVT. GIT (Green Infrastructure Typology) is suitable for use at the local scale with 50-100m data grids, but needs calibration, especially for mixed types of land cover. While GEDI's 30m diameter shots can help, the fundamental problem is that high-resolution GIT-type classification currently requires expensive airborne remote sensing data. Replacement by low-cost space-borne LIDAR data would be a game changer.  

Aside from making use of the perhaps unexpected high accuracy of GEDI-generated digital terrain models over water (PDF) and in urban areas (PDF), and recognising the average GEDI geolocation error of about 8 to 11 m (PDF), the main GEDI focus is not the shot level but the 1000m mesh level, primarily to estimate mean biomass using models. 

It is estimated that globally and certainly in the northern hemisphere, most of the Earth's 1000m grid cells will be intercepted by two GEDI beam swathes. There are numerous options for interpolating point data to continuous grids depending on the nature of the dataset. As a starting point, a kriging methodology similar to the one outlined in the GEDI Level 3 Algorithm Theoretical Basis Document (PDF) is being used. This is a good option for natural and continuous forest types but probably not for mixed land cover types (called fragmented or urban/suburban forests), vegetated built-up areas and nature-based infrastructure. These will probably need advanced interpolation methods such as co-kriging or fusion with auxiliary data having continuous spatial coverage such as satellite imagery and possibly OpenStreetMap data.

Microclimate classification is becoming the basis for infrastructure project design and selection. If and how GEDI addresses meshing will therefore probably impact the way this classification could be carried out. GEDI height data currently do not discriminate trees from buildings so the Global Forest Canopy Height service for example masks urban areas on the basis of a vegetation index and the gridded Global Human Settlement Layer (GHLS) built-up area dataset. Planned are more sophisticated approaches that would normally include a broader range of land cover masks (e.g., inland water, barren land, imperviousness, and wetland). Filtering data for built infrastructure and vegetation may be another approach. Similar methods are probably applicable for climate mitigation carbon accounting for project planning and certification and policy making at the non-forest project-level, where trees accompanying built infrastructure are an important biomass contributor (PDF).

A similar discussion can take place for project-level biodiversity accounting based on earth observation. So after making a carbon account for one of our vegetation-rich villages,  Jura Mountains mapping may turn its project-scale ecology accounting towards biodiversity accounting, perhaps in terms of functional diversity (PDF).


QGIS

Manipulating GEDI data using QGIS is fairly straightforward. As explained, the GEDI subsetter should be copied to the conda enviroment. It is then run using a terminal command of the form:

  • python GEDI_Subsetter.py --dir <insert input directory with GEDI files here> --roi <UpperLeftLatitude,UpperLeftLongitude,LowerRightLatitude,LowerRightLongitude>

which for our simple Windows10 environment translates to :

  • python GEDI_Subsetter.py --dir C:\Users\admin\ --roi 46.75,6.0,46.35,6.5

This creates a geojson .json file with the same name as the GEDI .h5 file in a directory call "output". The json file can be loaded into QGIS.

In QGIS, for GEDI Level 2B data, shots can be classified by biophysical  attributes such as Plant Area Index (pai), Total Canopy Cover (cover), and Foliage Height Diversity (fhd_normal) and elevation (elev_lowestmode). An archeological investigation provides an illustration. The QGIS clasification is done by choosing layer properties and setting the Symbology to "Categorized" and the "Value" to say rh100, selecting a "Color ramp", and clicking "Classify".

It is useful to have custom layers such as a Digital Surface Model of terrain available in QGIS.  So install the OpenLayers plugin in QGIS and for Windows, as explained, edit say the osm.py file

  • c:\Users\admin\AppData\Roaming\QGIS\QGIS3\profiles\default\python\plugin\openlayers_plugin\weblayers

replacing say the url for the Humanitarian Data Model with a custom url (something like http://map.peterboswell.net/openarea/{z}({x}/{y}.png in our case).

Editing files to set the transparency of the cusom layer in QGIS is a bit fiddly. Best to simply drag the GEDI layer above any custom or standard image layers such a Bing Satellite. 

An R package is also available for manipuating GEDI products.


Diversity

The functional diversity of forests may better predict ecosystem functioning than a species-level approach. Morphological traits can be used to describe a forest canopy's architecture in terms of different aspects of functional diversity, namely functional richness, divergence and evenness.

Our immediate interest is to see if there is a way to identify small-scale, localised forest disturbance using morphological traits estimated using GEDI data.

Starting with Level_2B data, we give on the forest map for the Series 3 transects, the laser shot points as RGB colour composites of three normalised GEDI morphological traits:

  • Total Canopy Cover (cover)
  • Total Plant Area Index (pai)
  • Foliage Height Diversity (fhd)

As for  the ternary mapping of morphological forest traits which used  the normalised:

  • Canopy Height (CH, canopy top to ground distance) 
  • PAI (projected plant area per horizontal ground area)
  • FHD

our map's MORPHOLOGY track is coloured as follows:

  • Red:    COVER <0.5  FHD <0.5  PAI>0.5 
  • Blue:   COVER <0.5  FHD >0.5  PAI<0.5
  • Green: COVER >0.3  FHD >0.5  PAI<0.5

The transects were measured in mid-winter with snow-on conditions so should be taken as a proof-of-concept at this stage.

It can be seen that grass-covered areas are blue,  and forested areas generally pink since vegetation is exposed. Leaf-on spring-to-summer transects are being developed. 

The composite RGB morphology track can be aligned to the accompanying PAVD plot by rotating and zooming the map appropriately.

The technique is straightforward:  

  1. Downloaded GEDI scans are converted to .geojson files using GEDI_Subsetter.py (see above). In QGIS, the .geojson is reprojected to EPSG3857 and split into two files (high quality and low quality) using the QGIS's Data Management Tools > Split Vector Layers to split according to the l2b_quality_flag (value  = 1 for high quality). 
  2. The high-quality .geojson is then rasterized into three separate files using Raster > Conversion > Rasterize using a horizontal resolution of 25 and a vertical resolution of 50, with burn value fields cover, fhd and pai.
  3. A virtual raster is then built from the three rasterized files using Raster > Miscellaneous > Build Virtual Raster  with "each input file placed in a separate  band" checked. 
  4. The merged file's layer properties are set by setting each band's min and max to 0 and 1 and adding a transparency  pixel. 
  5. The  merged RGB composite  is then exported as a .tif file and the QGIS "Generate XYZ tiles" tool is used to create tiles.

The current plan is to use both GEDI Level_2A and Level_2B data to generate composite RGB point tracks for the various GEDI structural traits that are available as GEDI  Level_2A and Level_2B products (see ATBD - PDF2A Levels - PDF; 2B Levels - PDF):

  • Total Canopy Cover or Canopy Cover Fraction (cover)
  • Cumulative canopy cover vertical profile or  Vertical Cover profile (cover-z)
  • Height above ground of the received waveform signal start (Level_2A: rh101; Level_2B: rh100)
  • Total Plant Area Index (pai
  • Plant Area Index profile (pai-z)
  • Plant Area Volume Density profile (pavd-z)
  • Foliage Height Diversity (fhd)

Other structural traits that are routinely investigated include:

  • Canopy Height 
  • Relative Heights (RH) as percentiles of the vertical distribution of canopy points at 25, 50, 75 and 95%
  • Canopy Ratio (canopy depth to canopy height given by [RH98 − RH25]/RH98 where RH refers to the Relative Height at the 25% and 98% percentiles)
  • PAI at 10m height intervals, notably PAI0-10, PAI10-20, etc 
  • Understory Density (PAI0-10)

When it comes to functional richness as opposed to say functional diversity, the following traits probably have the strongest and most independent contributions (Schneider et al, 2020):

  • Canopy Ratio ([RH98 − RH25]/RH98)
  • PAI (pai)
  • Understory Density (PAI0-10)

Terrain roughness is clearly an issue, especially topographic effects in steep terrain. For example, our RGB composite tracks overlay a high-resolution surface model (SwissSURFACE) and high resolution arial imagery (SwissIMAGE) to show that points that strike steep, forested areas, possibly with deep undersory vegetation, display unusually low PAI values. Any conclusion along these lines must of course recognise geolocation errors that arise not only for the GEDI tracks but also for the aerial imagery.

In passing, we note that the SwissIMAGE high-resolution imagery that is now publically available can be displayed not only in the JOSM OpenStreetMap map editor and other editors, but also in Leaflet and OpenLayers map layers using respectively: 

  • var swissimageLayer = L.tileLayer.wms('https://wms.geo.admin.ch/?FORMAT=image/vnd.jpeg-png8&TRANSPARENT=TRUE&VERSION=1.3.0&SERVICE=WMS&REQUEST=GetMap&LAYERS=ch.swisstopo.images-swissimage&STYLES=&CRS=EPSG:3857',{layers:'ch.swisstopo.images-swissimage',format:'image/png'});
  • var osmAerial = new ol.layer.Tile({source: new ol.source.TileWMS({crossOrigin: null, url: 'https://wms.geo.admin.ch/', params: {'LAYERS': 'ch.swisstopo.images-swissimage','FORMAT': 'image/jpeg',},serverType: 'mapserver',})});

April 2021


Email: peter@peterboswell.com
Phone: +41792989666
Skype: live:petergboswell
WhatsApp: PeterBoswellcom
Messenger: petergboswell
Telegram: PeterBoswellBot
LinkedIn
Facebook
Twitter
OneName