{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Getting Started with GEDI L2A Version 2 Data in Python\n", "### This tutorial demonstrates how to work with the Elevation and Height Metrics ([GEDI02_A.002](https://doi.org/10.5067/GEDI/GEDI02_A.002)) data product.\n", "The Global Ecosystem Dynamics Investigation ([GEDI](https://lpdaac.usgs.gov/data/get-started-data/collection-overview/missions/gedi-overview/)) mission aims to characterize ecosystem structure and dynamics to enable radically improved quantification and understanding of the Earth's carbon cycle and biodiversity. The GEDI instrument produces high resolution laser ranging observations of the 3-dimensional structure of the Earth. GEDI is attached to the International Space Station and collects data globally between 51.6 N and 51.6 S latitudes at the highest resolution and densest sampling of any light detection and ranging (lidar) instrument in orbit to date. The Land Processes Distributed Active Archive Center (LP DAAC) distributes the GEDI Level 1 and Level 2 Version 1 and Version 2 products. The L1B and L2 GEDI products are archived and distributed in the HDF-EOS5 file format. \n", "\n", "---\n", "## Use Case Example: \n", "This tutorial was developed using an example use case for a project being completed by the National Park Service. **The goal of the project is to use GEDI L2A Version 2 data to observe tree canopy height and profile over Redwood National Park in northern California.** \n", "\n", "This tutorial will show how to use Python to open GEDI L2A Version 2 files, visualize the sub-orbit of GEDI points (shots), subset to a region of interest, visualize GEDI canopy height, and export subsets of GEDI science dataset (SDS) layers as GeoJSON files that can be loaded into GIS and/or Remote Sensing software programs. \n", "\n", "- [Redwood National Park GeoJSON](https://github.com/nasa/GEDI-Data-Resources/data/RedwoodNP.geojson) \n", " - Contains the administrative boundary for Redwood National Park, available from: [Administrative Boundaries of National Park System Units 12/31/2017 - National Geospatial Data Asset (NGDA) NPS National Parks Dataset](https://irma.nps.gov/DataStore/DownloadFile/594958)\n", "- [waveform.csv](https://github.com/nasa/GEDI-Data-Resources/data/waveform.csv) \n", " \n", "*** \n", "### Data Used in the Example: \n", "- **GEDI L2A Elevation and Height Metrics Data Global Footprint Level - [GEDI02_A.002](https://doi.org/10.5067/GEDI/GEDI02_A.002)**\n", " - _The purpose of the L2A dataset is to provide waveform interpretation and extracted products from each GEDI waveform. This 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._\n", " - **Science Dataset (SDS) layers:**\n", " - /rh \n", " - /elev_lowestmode \n", " - /elev_highestreturn \n", " - /lat_lowestmode \n", " - /lon_lowestmode \n", " - /shot_number \n", " - /quality_flag \n", " - /digital_elevation_model \n", " - /degrade_flag \n", " - /sensitivity \n", " - /selected_algorithm\n", "*** \n", "# Topics Covered:\n", "1. [**Get Started**](#getstarted) \n", " 1.1 Import Packages \n", " 1.2 Set Up the Working Environment and Retrieve Files \n", "2. [**Import and Interpret Data**](#importinterpret) \n", " 2.1 Open a GEDI HDF5 File and Read File Metadata \n", " 2.2 Read SDS Metadata and Subset by Beam \n", "3. [**Visualize a GEDI Sub-Orbit**](#visualizeorbit) \n", " 3.1 Subset by Layer and Create a Geodataframe \n", " 3.2 Visualize a Geodataframe\n", "4. [**Work with GEDI L2A Data**](#L2A) \n", " 4.1 Import and Extract Specific Shots \n", " 4.2 Plot Relative Height Metrics \n", " 4.3 Combine RH Metrics and Waveforms \n", " 4.4 Select Data from non-Default Algorithm \n", "5. [**Plot Transects**](#plottransects) \n", " 5.1 Quality Filtering \n", " 5.2 Plot Beam Transects \n", " 5.3 Subset Beam Transects \n", " 5.4 Plot RH Metrics Transects \n", "6. [**Spatial Visualization**](#spatialvisualization) \n", " 6.1 Import, Subset, and Quality Filter all Beams \n", " 6.2 Spatial Subsetting \n", " 6.3 Visualize All Beams: Canopy Height and Elevation \n", "7. [**Export Subsets as GeoJSON Files**](#exportgeojson) \n", "***\n", "# Before Starting this Tutorial:\n", "\n", "This tutorial requires a compatible Python Environment and GEDI L2A Version 2 granule from June 19, 2019 (orbit 02932, sub-orbit `02`) to download. To setup the Python environment and download the file, follow the steps in sections 1 and 2 of the [set-up instruction](https://github.com/nasa/GEDI-Data-Resources/Setup/setup_instructions.md). \n", "\n", "## Source Code used to Generate this Tutorial:\n", "The repository containing all of the required files is located at: https://github.com/nasa/GEDI-Data-Resources \n", "\n", "
\n", "NOTE: This tutorial was developed for GEDI L2A Version 2 files and should only be used for those products.
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "# 1. Get Started " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1.1 Import Packages \n", "#### Import the required packages and set the input/working directory to run this Jupyter Notebook locally." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/javascript": "(function(root) {\n function now() {\n return new Date();\n }\n\n const force = true;\n const py_version = '3.7.2'.replace('rc', '-rc.').replace('.dev', '-dev.');\n const reloading = false;\n const Bokeh = root.Bokeh;\n\n // Set a timeout for this load but only if we are not already initializing\n if (typeof (root._bokeh_timeout) === \"undefined\" || (force || !root._bokeh_is_initializing)) {\n root._bokeh_timeout = Date.now() + 5000;\n root._bokeh_failed_load = false;\n }\n\n function run_callbacks() {\n try {\n root._bokeh_onload_callbacks.forEach(function(callback) {\n if (callback != null)\n callback();\n });\n } finally {\n delete root._bokeh_onload_callbacks;\n }\n console.debug(\"Bokeh: all callbacks have finished\");\n }\n\n function load_libs(css_urls, js_urls, js_modules, js_exports, callback) {\n if (css_urls == null) css_urls = [];\n if (js_urls == null) js_urls = [];\n if (js_modules == null) js_modules = [];\n if (js_exports == null) js_exports = {};\n\n root._bokeh_onload_callbacks.push(callback);\n\n if (root._bokeh_is_loading > 0) {\n // Don't load bokeh if it is still initializing\n console.debug(\"Bokeh: BokehJS is being loaded, scheduling callback at\", now());\n return null;\n } else if (js_urls.length === 0 && js_modules.length === 0 && Object.keys(js_exports).length === 0) {\n // There is nothing to load\n run_callbacks();\n return null;\n }\n\n function on_load() {\n root._bokeh_is_loading--;\n if (root._bokeh_is_loading === 0) {\n console.debug(\"Bokeh: all BokehJS libraries/stylesheets loaded\");\n run_callbacks()\n }\n }\n window._bokeh_on_load = on_load\n\n function on_error(e) {\n const src_el = e.srcElement\n console.error(\"failed to load \" + (src_el.href || src_el.src));\n }\n\n const skip = [];\n if (window.requirejs) {\n window.requirejs.config({'packages': {}, 'paths': {}, 'shim': {}});\n root._bokeh_is_loading = css_urls.length + 0;\n } else {\n root._bokeh_is_loading = css_urls.length + js_urls.length + js_modules.length + Object.keys(js_exports).length;\n }\n\n const existing_stylesheets = []\n const links = document.getElementsByTagName('link')\n for (let i = 0; i < links.length; i++) {\n const link = links[i]\n if (link.href != null) {\n existing_stylesheets.push(link.href)\n }\n }\n for (let i = 0; i < css_urls.length; i++) {\n const url = css_urls[i];\n const escaped = encodeURI(url)\n if (existing_stylesheets.indexOf(escaped) !== -1) {\n on_load()\n continue;\n }\n const element = document.createElement(\"link\");\n element.onload = on_load;\n element.onerror = on_error;\n element.rel = \"stylesheet\";\n element.type = \"text/css\";\n element.href = url;\n console.debug(\"Bokeh: injecting link tag for BokehJS stylesheet: \", url);\n document.body.appendChild(element);\n } var existing_scripts = []\n const scripts = document.getElementsByTagName('script')\n for (let i = 0; i < scripts.length; i++) {\n var script = scripts[i]\n if (script.src != null) {\n existing_scripts.push(script.src)\n }\n }\n for (let i = 0; i < js_urls.length; i++) {\n const url = js_urls[i];\n const escaped = encodeURI(url)\n if (skip.indexOf(escaped) !== -1 || existing_scripts.indexOf(escaped) !== -1) {\n if (!window.requirejs) {\n on_load();\n }\n continue;\n }\n const element = document.createElement('script');\n element.onload = on_load;\n element.onerror = on_error;\n element.async = false;\n element.src = url;\n console.debug(\"Bokeh: injecting script tag for BokehJS library: \", url);\n document.head.appendChild(element);\n }\n for (let i = 0; i < js_modules.length; i++) {\n const url = js_modules[i];\n const escaped = encodeURI(url)\n if (skip.indexOf(escaped) !== -1 || existing_scripts.indexOf(escaped) !== -1) {\n if (!window.requirejs) {\n on_load();\n }\n continue;\n }\n var element = document.createElement('script');\n element.onload = on_load;\n element.onerror = on_error;\n element.async = false;\n element.src = url;\n element.type = \"module\";\n console.debug(\"Bokeh: injecting script tag for BokehJS library: \", url);\n document.head.appendChild(element);\n }\n for (const name in js_exports) {\n const url = js_exports[name];\n const escaped = encodeURI(url)\n if (skip.indexOf(escaped) >= 0 || root[name] != null) {\n if (!window.requirejs) {\n on_load();\n }\n continue;\n }\n var element = document.createElement('script');\n element.onerror = on_error;\n element.async = false;\n element.type = \"module\";\n console.debug(\"Bokeh: injecting script tag for BokehJS library: \", url);\n element.textContent = `\n import ${name} from \"${url}\"\n window.${name} = ${name}\n window._bokeh_on_load()\n `\n document.head.appendChild(element);\n }\n if (!js_urls.length && !js_modules.length) {\n on_load()\n }\n };\n\n function inject_raw_css(css) {\n const element = document.createElement(\"style\");\n element.appendChild(document.createTextNode(css));\n document.body.appendChild(element);\n }\n\n const js_urls = [\"https://cdn.holoviz.org/panel/1.6.3/dist/bundled/reactiveesm/es-module-shims@^1.10.0/dist/es-module-shims.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-3.7.2.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-gl-3.7.2.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-widgets-3.7.2.min.js\", \"https://cdn.bokeh.org/bokeh/release/bokeh-tables-3.7.2.min.js\", \"https://cdn.holoviz.org/panel/1.6.3/dist/panel.min.js\", \"https://cdn.jsdelivr.net/npm/@holoviz/geoviews@1.14.0/dist/geoviews.min.js\"];\n const js_modules = [];\n const js_exports = {};\n const css_urls = [];\n const inline_js = [ function(Bokeh) {\n Bokeh.set_log_level(\"info\");\n },\nfunction(Bokeh) {} // ensure no trailing comma for IE\n ];\n\n function run_inline_js() {\n if ((root.Bokeh !== undefined) || (force === true)) {\n for (let i = 0; i < inline_js.length; i++) {\n try {\n inline_js[i].call(root, root.Bokeh);\n } catch(e) {\n if (!reloading) {\n throw e;\n }\n }\n }\n // Cache old bokeh versions\n if (Bokeh != undefined && !reloading) {\n var NewBokeh = root.Bokeh;\n if (Bokeh.versions === undefined) {\n Bokeh.versions = new Map();\n }\n if (NewBokeh.version !== Bokeh.version) {\n Bokeh.versions.set(NewBokeh.version, NewBokeh)\n }\n root.Bokeh = Bokeh;\n }\n } else if (Date.now() < root._bokeh_timeout) {\n setTimeout(run_inline_js, 100);\n } else if (!root._bokeh_failed_load) {\n console.log(\"Bokeh: BokehJS failed to load within specified timeout.\");\n root._bokeh_failed_load = true;\n }\n root._bokeh_is_initializing = false\n }\n\n function load_or_wait() {\n // Implement a backoff loop that tries to ensure we do not load multiple\n // versions of Bokeh and its dependencies at the same time.\n // In recent versions we use the root._bokeh_is_initializing flag\n // to determine whether there is an ongoing attempt to initialize\n // bokeh, however for backward compatibility we also try to ensure\n // that we do not start loading a newer (Panel>=1.0 and Bokeh>3) version\n // before older versions are fully initialized.\n if (root._bokeh_is_initializing && Date.now() > root._bokeh_timeout) {\n // If the timeout and bokeh was not successfully loaded we reset\n // everything and try loading again\n root._bokeh_timeout = Date.now() + 5000;\n root._bokeh_is_initializing = false;\n root._bokeh_onload_callbacks = undefined;\n root._bokeh_is_loading = 0\n console.log(\"Bokeh: BokehJS was loaded multiple times but one version failed to initialize.\");\n load_or_wait();\n } else if (root._bokeh_is_initializing || (typeof root._bokeh_is_initializing === \"undefined\" && root._bokeh_onload_callbacks !== undefined)) {\n setTimeout(load_or_wait, 100);\n } else {\n root._bokeh_is_initializing = true\n root._bokeh_onload_callbacks = []\n const bokeh_loaded = root.Bokeh != null && (root.Bokeh.version === py_version || (root.Bokeh.versions !== undefined && root.Bokeh.versions.has(py_version)));\n if (!reloading && !bokeh_loaded) {\n if (root.Bokeh) {\n root.Bokeh = undefined;\n }\n console.debug(\"Bokeh: BokehJS not loaded, scheduling load and callback at\", now());\n }\n load_libs(css_urls, js_urls, js_modules, js_exports, function() {\n console.debug(\"Bokeh: BokehJS plotting callback run at\", now());\n run_inline_js();\n });\n }\n }\n // Give older versions of the autoload script a head-start to ensure\n // they initialize before we start loading newer version.\n setTimeout(load_or_wait, 100)\n}(window));", "application/vnd.holoviews_load.v0+json": "" }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/javascript": "\nif ((window.PyViz === undefined) || (window.PyViz instanceof HTMLElement)) {\n window.PyViz = {comms: {}, comm_status:{}, kernels:{}, receivers: {}, plot_index: []}\n}\n\n\n function JupyterCommManager() {\n }\n\n JupyterCommManager.prototype.register_target = function(plot_id, comm_id, msg_handler) {\n if (window.comm_manager || ((window.Jupyter !== undefined) && (Jupyter.notebook.kernel != null))) {\n var comm_manager = window.comm_manager || Jupyter.notebook.kernel.comm_manager;\n comm_manager.register_target(comm_id, function(comm) {\n comm.on_msg(msg_handler);\n });\n } else if ((plot_id in window.PyViz.kernels) && (window.PyViz.kernels[plot_id])) {\n window.PyViz.kernels[plot_id].registerCommTarget(comm_id, function(comm) {\n comm.onMsg = msg_handler;\n });\n } else if (typeof google != 'undefined' && google.colab.kernel != null) {\n google.colab.kernel.comms.registerTarget(comm_id, (comm) => {\n var messages = comm.messages[Symbol.asyncIterator]();\n function processIteratorResult(result) {\n var message = result.value;\n console.log(message)\n var content = {data: message.data, comm_id};\n var buffers = []\n for (var buffer of message.buffers || []) {\n buffers.push(new DataView(buffer))\n }\n var metadata = message.metadata || {};\n var msg = {content, buffers, metadata}\n msg_handler(msg);\n return messages.next().then(processIteratorResult);\n }\n return messages.next().then(processIteratorResult);\n })\n }\n }\n\n JupyterCommManager.prototype.get_client_comm = function(plot_id, comm_id, msg_handler) {\n if (comm_id in window.PyViz.comms) {\n return window.PyViz.comms[comm_id];\n } else if (window.comm_manager || ((window.Jupyter !== undefined) && (Jupyter.notebook.kernel != null))) {\n var comm_manager = window.comm_manager || Jupyter.notebook.kernel.comm_manager;\n var comm = comm_manager.new_comm(comm_id, {}, {}, {}, comm_id);\n if (msg_handler) {\n comm.on_msg(msg_handler);\n }\n } else if ((plot_id in window.PyViz.kernels) && (window.PyViz.kernels[plot_id])) {\n var comm = window.PyViz.kernels[plot_id].connectToComm(comm_id);\n comm.open();\n if (msg_handler) {\n comm.onMsg = msg_handler;\n }\n } else if (typeof google != 'undefined' && google.colab.kernel != null) {\n var comm_promise = google.colab.kernel.comms.open(comm_id)\n comm_promise.then((comm) => {\n window.PyViz.comms[comm_id] = comm;\n if (msg_handler) {\n var messages = comm.messages[Symbol.asyncIterator]();\n function processIteratorResult(result) {\n var message = result.value;\n var content = {data: message.data};\n var metadata = message.metadata || {comm_id};\n var msg = {content, metadata}\n msg_handler(msg);\n return messages.next().then(processIteratorResult);\n }\n return messages.next().then(processIteratorResult);\n }\n })\n var sendClosure = (data, metadata, buffers, disposeOnDone) => {\n return comm_promise.then((comm) => {\n comm.send(data, metadata, buffers, disposeOnDone);\n });\n };\n var comm = {\n send: sendClosure\n };\n }\n window.PyViz.comms[comm_id] = comm;\n return comm;\n }\n window.PyViz.comm_manager = new JupyterCommManager();\n \n\n\nvar JS_MIME_TYPE = 'application/javascript';\nvar HTML_MIME_TYPE = 'text/html';\nvar EXEC_MIME_TYPE = 'application/vnd.holoviews_exec.v0+json';\nvar CLASS_NAME = 'output';\n\n/**\n * Render data to the DOM node\n */\nfunction render(props, node) {\n var div = document.createElement(\"div\");\n var script = document.createElement(\"script\");\n node.appendChild(div);\n node.appendChild(script);\n}\n\n/**\n * Handle when a new output is added\n */\nfunction handle_add_output(event, handle) {\n var output_area = handle.output_area;\n var output = handle.output;\n if ((output.data == undefined) || (!output.data.hasOwnProperty(EXEC_MIME_TYPE))) {\n return\n }\n var id = output.metadata[EXEC_MIME_TYPE][\"id\"];\n var toinsert = output_area.element.find(\".\" + CLASS_NAME.split(' ')[0]);\n if (id !== undefined) {\n var nchildren = toinsert.length;\n var html_node = toinsert[nchildren-1].children[0];\n html_node.innerHTML = output.data[HTML_MIME_TYPE];\n var scripts = [];\n var nodelist = html_node.querySelectorAll(\"script\");\n for (var i in nodelist) {\n if (nodelist.hasOwnProperty(i)) {\n scripts.push(nodelist[i])\n }\n }\n\n scripts.forEach( function (oldScript) {\n var newScript = document.createElement(\"script\");\n var attrs = [];\n var nodemap = oldScript.attributes;\n for (var j in nodemap) {\n if (nodemap.hasOwnProperty(j)) {\n attrs.push(nodemap[j])\n }\n }\n attrs.forEach(function(attr) { newScript.setAttribute(attr.name, attr.value) });\n newScript.appendChild(document.createTextNode(oldScript.innerHTML));\n oldScript.parentNode.replaceChild(newScript, oldScript);\n });\n if (JS_MIME_TYPE in output.data) {\n toinsert[nchildren-1].children[1].textContent = output.data[JS_MIME_TYPE];\n }\n output_area._hv_plot_id = id;\n if ((window.Bokeh !== undefined) && (id in Bokeh.index)) {\n window.PyViz.plot_index[id] = Bokeh.index[id];\n } else {\n window.PyViz.plot_index[id] = null;\n }\n } else if (output.metadata[EXEC_MIME_TYPE][\"server_id\"] !== undefined) {\n var bk_div = document.createElement(\"div\");\n bk_div.innerHTML = output.data[HTML_MIME_TYPE];\n var script_attrs = bk_div.children[0].attributes;\n for (var i = 0; i < script_attrs.length; i++) {\n toinsert[toinsert.length - 1].childNodes[1].setAttribute(script_attrs[i].name, script_attrs[i].value);\n }\n // store reference to server id on output_area\n output_area._bokeh_server_id = output.metadata[EXEC_MIME_TYPE][\"server_id\"];\n }\n}\n\n/**\n * Handle when an output is cleared or removed\n */\nfunction handle_clear_output(event, handle) {\n var id = handle.cell.output_area._hv_plot_id;\n var server_id = handle.cell.output_area._bokeh_server_id;\n if (((id === undefined) || !(id in PyViz.plot_index)) && (server_id !== undefined)) { return; }\n var comm = window.PyViz.comm_manager.get_client_comm(\"hv-extension-comm\", \"hv-extension-comm\", function () {});\n if (server_id !== null) {\n comm.send({event_type: 'server_delete', 'id': server_id});\n return;\n } else if (comm !== null) {\n comm.send({event_type: 'delete', 'id': id});\n }\n delete PyViz.plot_index[id];\n if ((window.Bokeh !== undefined) & (id in window.Bokeh.index)) {\n var doc = window.Bokeh.index[id].model.document\n doc.clear();\n const i = window.Bokeh.documents.indexOf(doc);\n if (i > -1) {\n window.Bokeh.documents.splice(i, 1);\n }\n }\n}\n\n/**\n * Handle kernel restart event\n */\nfunction handle_kernel_cleanup(event, handle) {\n delete PyViz.comms[\"hv-extension-comm\"];\n window.PyViz.plot_index = {}\n}\n\n/**\n * Handle update_display_data messages\n */\nfunction handle_update_output(event, handle) {\n handle_clear_output(event, {cell: {output_area: handle.output_area}})\n handle_add_output(event, handle)\n}\n\nfunction register_renderer(events, OutputArea) {\n function append_mime(data, metadata, element) {\n // create a DOM node to render to\n var toinsert = this.create_output_subarea(\n metadata,\n CLASS_NAME,\n EXEC_MIME_TYPE\n );\n this.keyboard_manager.register_events(toinsert);\n // Render to node\n var props = {data: data, metadata: metadata[EXEC_MIME_TYPE]};\n render(props, toinsert[0]);\n element.append(toinsert);\n return toinsert\n }\n\n events.on('output_added.OutputArea', handle_add_output);\n events.on('output_updated.OutputArea', handle_update_output);\n events.on('clear_output.CodeCell', handle_clear_output);\n events.on('delete.Cell', handle_clear_output);\n events.on('kernel_ready.Kernel', handle_kernel_cleanup);\n\n OutputArea.prototype.register_mime_type(EXEC_MIME_TYPE, append_mime, {\n safe: true,\n index: 0\n });\n}\n\nif (window.Jupyter !== undefined) {\n try {\n var events = require('base/js/events');\n var OutputArea = require('notebook/js/outputarea').OutputArea;\n if (OutputArea.prototype.mime_types().indexOf(EXEC_MIME_TYPE) == -1) {\n register_renderer(events, OutputArea);\n }\n } catch(err) {\n }\n}\n", "application/vnd.holoviews_load.v0+json": "" }, "metadata": {}, "output_type": "display_data" }, { "data": { "application/vnd.holoviews_exec.v0+json": "", "text/html": [ "
\n", "
\n", "
\n", "" ] }, "metadata": { "application/vnd.holoviews_exec.v0+json": { "id": "8d518a3b-6f6a-4a2b-a848-d9b01df9d2d0" } }, "output_type": "display_data" }, { "data": { "text/html": [ "\n", "
\n", "\n", "\n", "\n", " \n", " \n", "\n", "\n", " \n", " \n", "\n", "\n", "
\n" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import os\n", "import h5py\n", "import numpy as np\n", "import pandas as pd\n", "import geopandas as gp\n", "from shapely.geometry import Point\n", "import geoviews as gv\n", "from geoviews import opts, tile_sources as gvts\n", "import holoviews as hv\n", "gv.extension('bokeh', 'matplotlib')\n", "import shapely\n", "import warnings\n", "from shapely.errors import ShapelyDeprecationWarning\n", "warnings.filterwarnings(\"ignore\", category=ShapelyDeprecationWarning)\n", "\n", "os.chdir('../../') " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1.2 Set Up the Working Environment and Retrieve Files\n", "#### The input directory is defined as the current working directory. Note that you will need to have the Jupyter Notebook and example data (.h5 and .geojson) stored in this directory in order to execute the tutorial successfully." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "inDir = os.getcwd()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### You will need to download the file in order to execute this tutorial. Make sure to download the file into the `data` directory defined above.\n", "### Direct Link to file:\n", " - https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/GEDI02_A.002/GEDI02_A_2019170155833_O02932_02_T02267_02_003_01_V002/GEDI02_A_2019170155833_O02932_02_T02267_02_003_01_V002.h5 \n", "\n", "Alternatively, you can use `earthaccess` package to download the data. " ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "cf04af35f2e54387ada90d3fddb8a7ba", "version_major": 2, "version_minor": 0 }, "text/plain": [ "QUEUEING TASKS | : 0%| | 0/1 [00:00\n", "#### In this section, begin working with the [GEDI Level 2A Elevation and Height Metrics Data](https://doi.org/10.5067/GEDI/GEDI02_A.002). Begin by checking out the example L2A file metadata." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2.1 Open a GEDI HDF5 File and Read File Metadata \n", "#### Read the file using `h5py`. " ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'data/GEDI02_A_2019170155833_O02932_02_T02267_02_003_01_V002.h5'" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "L2A = f'data/{gediFiles[0]}'\n", "L2A" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### The standard format for GEDI filenames is as follows:\n", "> **GEDI02_A**: Product Short Name \n", "**2019170155833**: Julian Date and Time of Acquisition (YYYYDDDHHMMSS) \n", "**O02932**: Orbit Number \n", "**02**: Sub-Orbit Granule Number (1-4) \n", "**T02267**: Track Number (Reference Ground Track) \n", "**02**: Positioning and Pointing Determination System (PPDS) type (00 is predict, 01 rapid, 02 and higher is final) \n", "**003**: PGE Version Number \n", "**01**: Granule Production Version \n", "**V002**: Product Version " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Read in a GEDI HDF5 file using the `h5py` package." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "gediL2A = h5py.File(L2A, 'r') # Read file using h5py" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Navigate the HDF5 file below. " ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "['BEAM0000',\n", " 'BEAM0001',\n", " 'BEAM0010',\n", " 'BEAM0011',\n", " 'BEAM0101',\n", " 'BEAM0110',\n", " 'BEAM1000',\n", " 'BEAM1011',\n", " 'METADATA']" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "list(gediL2A.keys())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### The GEDI HDF5 file contains groups in which data and metadata are stored.\n", "#### First, the `METADATA` group contains the file-level metadata." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "['DatasetIdentification']" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "list(gediL2A['METADATA'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This contains useful information such as the creation date, PGEVersion, and VersionID. Below, print the file-level metadata attributes." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "PGEVersion\n", "VersionID\n", "abstract\n", "characterSet\n", "creationDate\n", "credit\n", "fileName\n", "language\n", "originatorOrganizationName\n", "purpose\n", "shortName\n", "spatialRepresentationType\n", "status\n", "topicCategory\n", "uuid\n" ] } ], "source": [ "for g in gediL2A['METADATA']['DatasetIdentification'].attrs: print(g) " ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The purpose of the L2A dataset is to provide waveform interpretation and extracted products from each GEDI waveform. This 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.\n" ] } ], "source": [ "print(gediL2A['METADATA']['DatasetIdentification'].attrs['purpose'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2.2 Read SDS Metadata and Subset by Beam " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### The GEDI instrument consists of 3 lasers producing a total of 8 beam ground transects. The eight remaining groups contain data for each of the eight GEDI beam transects. For additional information, be sure to check out: https://gedi.umd.edu/instrument/specifications/." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "['BEAM0000',\n", " 'BEAM0001',\n", " 'BEAM0010',\n", " 'BEAM0011',\n", " 'BEAM0101',\n", " 'BEAM0110',\n", " 'BEAM1000',\n", " 'BEAM1011']" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "beamNames = [g for g in gediL2A.keys() if g.startswith('BEAM')]\n", "beamNames" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### One useful piece of metadata to retrieve from each beam transect is whether it is a full power beam or a coverage beam. " ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "description\n" ] } ], "source": [ "for g in gediL2A['BEAM0000'].attrs: print(g)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "BEAM0000 is a Coverage beam\n", "BEAM0001 is a Coverage beam\n", "BEAM0010 is a Coverage beam\n", "BEAM0011 is a Coverage beam\n", "BEAM0101 is a Full power beam\n", "BEAM0110 is a Full power beam\n", "BEAM1000 is a Full power beam\n", "BEAM1011 is a Full power beam\n" ] } ], "source": [ "for b in beamNames: \n", " print(f\"{b} is a {gediL2A[b].attrs['description']}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Below, pick one of the full power beams that will be used to retrieve GEDI L2A relative height metrics in Section 3. " ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "beamNames = ['BEAM0110']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Identify all the objects in the GEDI HDF5 file below. \n", "Note: This step may take a while to complete." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "['BEAM0110/ancillary/l2a_alg_count',\n", " 'BEAM0110/beam',\n", " 'BEAM0110/channel',\n", " 'BEAM0110/degrade_flag',\n", " 'BEAM0110/delta_time',\n", " 'BEAM0110/digital_elevation_model',\n", " 'BEAM0110/digital_elevation_model_srtm',\n", " 'BEAM0110/elev_highestreturn',\n", " 'BEAM0110/elev_lowestmode',\n", " 'BEAM0110/elevation_bias_flag']" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "gediL2A_objs = []\n", "gediL2A.visit(gediL2A_objs.append) # Retrieve list of datasets\n", "gediSDS = [o for o in gediL2A_objs if isinstance(gediL2A[o], h5py.Dataset)] # Search for relevant SDS inside data file\n", "[i for i in gediSDS if beamNames[0] in i][0:10] # Print the first 10 datasets for selected beam" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "# 3. Visualize a GEDI Orbit \n", "#### In the section below, import GEDI L2A SDS layers into a `GeoPandas` GeoDataFrame for the beam specified above. \n", "#### Use the `lat_lowestmode` and `lon_lowestmode` to create a `shapely` point for each GEDI shot location. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3.1 Subset by Layer and Create a Geodataframe " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Read in the SDS and take a representative sample (every 100th shot) and append to lists, then use the lists to generate a `pandas` dataframe." ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
BeamShot NumberLongitudeLatitudeQuality Flag
0BEAM011029320600200419869-142.75569226.9238960
1BEAM011029320600200419969-142.73656726.9432420
2BEAM011029320600200420069-142.71743326.9625690
3BEAM011029320600200420169-142.69829426.9818910
4BEAM011029320600200420269-142.67913627.0011870
..................
1066BEAM011029320600200526469-80.19845251.7968580
1067BEAM011029320600200526569-80.11467751.7970290
1068BEAM011029320600200526669-80.03217351.7971730
1069BEAM011029320600200526769-79.94857151.7972450
1070BEAM011029320600200526869-79.86569851.7972460
\n", "

1071 rows × 5 columns

\n", "
" ], "text/plain": [ " Beam Shot Number Longitude Latitude Quality Flag\n", "0 BEAM0110 29320600200419869 -142.755692 26.923896 0\n", "1 BEAM0110 29320600200419969 -142.736567 26.943242 0\n", "2 BEAM0110 29320600200420069 -142.717433 26.962569 0\n", "3 BEAM0110 29320600200420169 -142.698294 26.981891 0\n", "4 BEAM0110 29320600200420269 -142.679136 27.001187 0\n", "... ... ... ... ... ...\n", "1066 BEAM0110 29320600200526469 -80.198452 51.796858 0\n", "1067 BEAM0110 29320600200526569 -80.114677 51.797029 0\n", "1068 BEAM0110 29320600200526669 -80.032173 51.797173 0\n", "1069 BEAM0110 29320600200526769 -79.948571 51.797245 0\n", "1070 BEAM0110 29320600200526869 -79.865698 51.797246 0\n", "\n", "[1071 rows x 5 columns]" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "lonSample, latSample, shotSample, qualitySample, beamSample = [], [], [], [], [] # Set up lists to store data\n", "\n", "# Open the SDS\n", "lats = gediL2A[f'{beamNames[0]}/lat_lowestmode'][()]\n", "lons = gediL2A[f'{beamNames[0]}/lon_lowestmode'][()]\n", "shots = gediL2A[f'{beamNames[0]}/shot_number'][()]\n", "quality = gediL2A[f'{beamNames[0]}/quality_flag'][()]\n", "\n", "# Take every 100th shot and append to list\n", "for i in range(len(shots)):\n", " if i % 100 == 0:\n", " shotSample.append(str(shots[i]))\n", " lonSample.append(lons[i])\n", " latSample.append(lats[i])\n", " qualitySample.append(quality[i])\n", " beamSample.append(beamNames[0])\n", " \n", "# Write all of the sample shots to a dataframe\n", "latslons = pd.DataFrame({'Beam': beamSample, 'Shot Number': shotSample, 'Longitude': lonSample, 'Latitude': latSample,\n", " 'Quality Flag': qualitySample})\n", "latslons" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Above is a dataframe containing columns describing the beam, shot number, lat/lon location, and quality information about each shot." ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "# Clean up variables that will no longer be needed\n", "del beamSample, quality, qualitySample, gediL2A_objs, latSample, lats, lonSample, lons, shotSample, shots " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Below, create an additional column called 'geometry' that contains a `shapely` point generated from each lat/lon location from the shot. " ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "source": [ "# Take the lat/lon dataframe and convert each lat/lon to a shapely point\n", "latslons['geometry'] = latslons.apply(lambda row: Point(row.Longitude, row.Latitude), axis=1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Next, convert to a `Geopandas` GeoDataFrame." ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0 POINT (-142.75569 26.9239)\n", "1 POINT (-142.73657 26.94324)\n", "2 POINT (-142.71743 26.96257)\n", "3 POINT (-142.69829 26.98189)\n", "4 POINT (-142.67914 27.00119)\n", " ... \n", "1066 POINT (-80.19845 51.79686)\n", "1067 POINT (-80.11468 51.79703)\n", "1068 POINT (-80.03217 51.79717)\n", "1069 POINT (-79.94857 51.79725)\n", "1070 POINT (-79.8657 51.79725)\n", "Name: geometry, Length: 1071, dtype: geometry" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Convert to a Geodataframe\n", "latslons = gp.GeoDataFrame(latslons)\n", "latslons = latslons.drop(columns=['Latitude','Longitude'])\n", "latslons['geometry']" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Pull out and plot an example `shapely` point below." ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "scrolled": true }, "outputs": [ { "data": { "image/svg+xml": [ "" ], "text/plain": [ "" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "latslons['geometry'][0]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3.2 Visualize a GeoDataFrame \n", "#### In this section, use the GeoDataFrame and the `geoviews` python package to spatially visualize the location of the GEDI shots on a basemap and import a GeoJSON file of the spatial region of interest for the use case example: Redwood National Park." ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [], "source": [ "# Define a function for visualizing GEDI points\n", "def pointVisual(features, vdims):\n", " return (gvts.EsriImagery * gv.Points(features, vdims=vdims).options(tools=['hover'], height=500, width=900, size=5, \n", " color='yellow', fontsize={'xticks': 10, 'yticks': 10, \n", " 'xlabel':16, 'ylabel': 16}))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Import a GeoJSON of Redwood National Park as an additional GeoDataFrame. Note that you will need to have downloaded the [GeoJSON](https://github.com/nasa/GEDI-Data-Resources/tree/main/data/RedwoodNP.geojson) from the bitbucket repo containing this tutorial and have it saved in the same directory as this Jupyter Notebook." ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [], "source": [ "redwoodNP = gp.GeoDataFrame.from_file('data/RedwoodNP.geojson') # Import GeoJSON as GeoDataFrame" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
GIS_LOC_IDUNIT_CODEGROUP_CODEUNIT_NAMEUNIT_TYPEMETA_MIDFLANDS_CODEDATE_EDITGIS_NOTESgeometry
0NoneREDWNoneRedwoodNational ParkNoneNoneNoneShifted 0.06 milesMULTIPOLYGON (((-124.01829 41.44539, -124.0184...
\n", "
" ], "text/plain": [ " GIS_LOC_ID UNIT_CODE GROUP_CODE UNIT_NAME UNIT_TYPE META_MIDF \\\n", "0 None REDW None Redwood National Park None \n", "\n", " LANDS_CODE DATE_EDIT GIS_NOTES \\\n", "0 None None Shifted 0.06 miles \n", "\n", " geometry \n", "0 MULTIPOLYGON (((-124.01829 41.44539, -124.0184... " ] }, "execution_count": 23, "metadata": {}, "output_type": "execute_result" } ], "source": [ "redwoodNP" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "data": { "image/svg+xml": [ "" ], "text/plain": [ "" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "redwoodNP['geometry'][0] # Plot GeoDataFrame" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Defining the vdims below will allow you to hover over specific shots and view information about them." ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "['Beam', 'Shot Number', 'Quality Flag']" ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Create a list of geodataframe columns to be included as attributes in the output map\n", "vdims = []\n", "for f in latslons:\n", " if f not in ['geometry']:\n", " vdims.append(f)\n", "vdims" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Below, combine a plot of the Redwood National Park Boundary (combine two `geoviews` plots using `*`) with the point visual mapping function defined above in order to plot (1) the representative GEDI shots, (2) the region of interest, and (3) a basemap layer. " ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "data": {}, "metadata": {}, "output_type": "display_data" }, { "data": { "application/vnd.holoviews_exec.v0+json": "", "text/html": [ "
\n", "
\n", "
\n", "" ], "text/plain": [ ":Overlay\n", " .Polygons.I :Polygons [Longitude,Latitude]\n", " .WMTS.I :WMTS [Longitude,Latitude]\n", " .Points.I :Points [Longitude,Latitude] (Beam,Shot Number,Quality Flag)" ] }, "execution_count": 26, "metadata": { "application/vnd.holoviews_exec.v0+json": { "id": "5e5e766d-a8fb-4764-b709-93761fe7c721" } }, "output_type": "execute_result" } ], "source": [ "# Call the function for plotting the GEDI points\n", "gv.Polygons(redwoodNP['geometry']).opts(line_color='red', color=None) * pointVisual(latslons, vdims = vdims)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Above is a good illustration of the new GEDI _Version 2_ sub-orbit granules (remember that GEDI _Version 1_ files are stored as one ISS orbit). One of the benefits of using geoviews is the interactive nature of the output plots. Use the tools to the right of the map above to zoom in and find the shots intersecting Redwood National Park. \n", "> (**HINT**: find where the orbit intersects the west coast of the United States)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Below is a screenshot of the region of interest:\n", "![GEDI_L2A_V2_Tutorial_1.png](../../img/GEDI_L2A_V2_Tutorial_1.png \"Sample of GEDI L1B shots in yellow (orbit 2932) plotted over Redwood National Park, USA.\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Side Note: Wondering what the 0's and 1's for `quality_flag` mean?" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Quality Flag: Flag simpilfying selection of most useful data\n" ] } ], "source": [ "print(f\"Quality Flag: {gediL2A[b]['quality_flag'].attrs['description']}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Above, 0 is poor quality and a quality_flag value of 1 indicates the laser shot meets criteria based on energy, sensitivity, amplitude, and real-time surface tracking quality. We will show an example of how to quality filter GEDI data in section 6.1.\n", "#### After finding one of the shots within Redwood NP, find the index for that shot number so that we can find the correct RH metrics to visualize in Section 4. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Each GEDI shot has a unique shot identifier (shot number) that is available within each data group of the product. The shot number is important to retain in any data subsetting as it will allow the user to link any shot record back to the original orbit data, and to link any shot and its data between the L1 and L2 products. The standard format for GEDI Shots is as follows:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Shot: 29320600200465601\n", "> **2932**: Orbit Number \n", "**06**: Beam Number \n", "**0**: Reserved for future use \n", "**02**: Sub-orbit Granule Number \n", "**004**: Minor frame number \n", "**65601**: Shot index " ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [], "source": [ "del latslons # No longer need the geodataframe used to visualize the full GEDI orbit" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "# 4. Work with GEDI L2A Data\n", "### Interpretation of RH Metrics\n", "#### The GEDI L2A data product provides relative height (RH) metrics, which are “lidar perceived” metrics that have the following characteristics:\n", " 1. RH100 = elev_highestreturn - elev_lowestmode\n", " 2. The RH metrics are intended for vegetated surfaces. Results over bare/water surfaces are still valid but may present some confusing results.\n", " 3. The lower RH metrics (e.g., RH10) will often have negative values, particularly in low canopy cover conditions. This is because a relatively high fraction of the waveform energy is from the ground and below elev_lowestmode. For example, if the ground return contains 30% of the energy, then RH1 through 15 are likely to be below 0 since half of the ground energy from the ground return is below the center of the ground return, which is used to determine the mean ground elevation in the footprint (elev_lowestmode).\n", "\n", "#### In this section, import and extract a specific GEDI L2A shot. From there, plot relative height (RH) metrics, and then combine with the L1B full waveform data. Note that we will continue to use `BEAM0110` from Section 3. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 4.1 Import and Extract Specific Shots\n", "#### Notice that there are thousands of datasets available in the GEDI L2A product. In the code blocks below, you will subset to just a few of the datasets available." ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "4320" ] }, "execution_count": 29, "metadata": {}, "output_type": "execute_result" } ], "source": [ "len(gediSDS)" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "['BEAM0110']" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" } ], "source": [ "beamNames" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Begin by subsetting just to the selected full power beam:" ] }, { "cell_type": "code", "execution_count": 31, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "540" ] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" } ], "source": [ "beamSDS = [g for g in gediSDS if beamNames[0] in g] # Subset to a single beam\n", "len(beamSDS)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### We will set the shot index used as an example from the [GEDI L1B Tutorial](https://github.com/nasa/GEDI-Data-Resources/blob/main/python/tutorials/GEDI_L1B_V2_Tutorial.ipynb) to show how to subset a single shot of GEDI L2A data." ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [], "source": [ "shot = 29320600200465601" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "np.int64(45732)" ] }, "execution_count": 33, "metadata": {}, "output_type": "execute_result" } ], "source": [ "index = np.where(gediL2A[f'{beamNames[0]}/shot_number'][()]==shot)[0][0] # Set the index for the shot identified above\n", "index" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 4.2 Plot Relative Height Metrics\n", "#### In section 4.2, import the relative height metrics and begin exploring how to plot them." ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [], "source": [ "rh = gediL2A[[g for g in beamSDS if g.endswith('/rh')][0]] # Relative Height Metrics" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Print the description for the `rh` dataset to better understand relative height (RH) metrics." ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "rh is Relative height metrics at 1 % interval\n" ] } ], "source": [ "print(f\"rh is {rh.attrs['description']}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### New with Version 2: \n", "#### The GEDI L2A product provides waveform processing results for multiple algorithm settings. Results for a default algorithm selection are provided in the root directory of the data product for each beam. In GEDI Version 1, the default for all shots was algorithm setting group 1. \n", "### In Version 2, the default algorithm is set to the best performing algorithm based on a number of characteristics for that specific shot. For information on selecting the non-default algorithm, see section 4.4 below. Below, bring in the `selected_algorithm` dataset to determine which algorithm (1-6) was set to the default for our shot. " ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "selected_algorithm is ID of algorithm selected as identifying the lowest non-noise mode\n" ] } ], "source": [ "algo = gediL2A[f'{beamNames[0]}/selected_algorithm'] # selected algorithm\n", "print(f\"selected_algorithm is {algo.attrs['description']}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Next, bring in other useful L2A datasets such as `lat_lowestmode` and `lon_lowestmode`." ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [], "source": [ "# Bring in the desired SDS\n", "lats = gediL2A[f'{beamNames[0]}/lat_lowestmode'] # Latitude\n", "lons = gediL2A[f'{beamNames[0]}/lon_lowestmode'] # Longitude" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Grab the location, relative height metrics, and algorithm for the shot defined above:" ] }, { "cell_type": "code", "execution_count": 38, "metadata": {}, "outputs": [], "source": [ "rhLat = lats[index]\n", "rhLon = lons[index]\n", "rhShot1 = rh[index]\n", "algoShot1 = algo[index]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Put everything together to identify the shot that we want to extract:" ] }, { "cell_type": "code", "execution_count": 39, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The shot is located at: 41.28541681552604, -124.03006394248702 (shot ID: 29320600200465601, index 45732) and is from beam BEAM0110.\n", "The selected algorithm is Algorithm Setting Group 1.\n" ] } ], "source": [ "print(f\"The shot is located at: {str(rhLat)}, {str(rhLon)} (shot ID: {shot}, index {index}) and is from beam {beamNames[0]}.\")\n", "print(f\"The selected algorithm is Algorithm Setting Group {str(algoShot1)}.\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### In order to plot the RH metrics, you also need to import the elevation recorded at `elev_lowestmode` (the first elevation recorded) and `elev_highestreturn` or the last elevation recorded for that waveform." ] }, { "cell_type": "code", "execution_count": 40, "metadata": {}, "outputs": [], "source": [ "# Grab the elevation recorded at the start and end of the RH metrics\n", "zElevation = gediL2A[[g for g in beamSDS if g.endswith('/elev_lowestmode')][0]][index] # Elevation\n", "zTop = gediL2A[[g for g in beamSDS if g.endswith('/elev_highestreturn')][0]][index] # Elevation Highest Return" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Below, convert canopy height to canopy elevation by adding the ground elevation to each RH 1% interval." ] }, { "cell_type": "code", "execution_count": 41, "metadata": {}, "outputs": [], "source": [ "rhShot = [z + zElevation for z in rhShot1] # To convert canopy height to canopy elevation, add the elevation to each value\n", "rh25 = rhShot[24] # 25% \n", "rh50 = rhShot[49] # 50% \n", "rh75 = rhShot[74] # 75% " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Next, use the `holoviews` (hv) package to make a line plot showing the elevation at each percentile for the selected shot." ] }, { "cell_type": "code", "execution_count": 42, "metadata": {}, "outputs": [ { "data": {}, "metadata": {}, "output_type": "display_data" }, { "data": { "application/vnd.holoviews_exec.v0+json": "", "text/html": [ "
\n", "
\n", "
\n", "" ], "text/plain": [ ":Curve [x] (y)" ] }, "execution_count": 42, "metadata": { "application/vnd.holoviews_exec.v0+json": { "id": "3c90a0b3-bb76-4a84-a7d1-e2a3e1e73db0" } }, "output_type": "execute_result" } ], "source": [ "rhVis = hv.Curve(rhShot, label=f'Selected Algorithm (a{str(algoShot1)})')\n", "rhVis = rhVis.opts(color='black', tools=['hover'], height=500, width=400, title='GEDI L2A Relative Height Metrics', \n", " xlabel='Percent Energy Returned', ylabel='Elevation (m)', xlim=(0,100),ylim=(np.min(rhShot),np.max(rhShot)), \n", " fontsize={'title':14, 'xlabel':16, 'ylabel': 16, 'legend': 14, 'xticks':12, 'yticks':12}, line_width=3.5)\n", "rhVis" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Congratulations! You have plotted your first relative cumulative RH profile.\n", "#### The `rh` dataset stores the relative height metrics at 1% intervals, and thus each shot contains 101 values representing rh at 0-100%." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Next, add some additional context to the cumulative RH profile, including metrics such as the `elev_lowestmode` (ground elevation), the `elev_highestreturn` (highest reflecting surface height), and the relative height metrics at each quartile." ] }, { "cell_type": "code", "execution_count": 43, "metadata": {}, "outputs": [], "source": [ "# Create plots for L2A Metrics\n", "zX = [0,100] # set up list from 0 to 100 to create the line\n", "zY = [zElevation, zElevation] # ground elevation\n", "zT = [zTop, zTop] # highest return\n", "\n", "# Set up plots for each of the desired values\n", "zVis = hv.Curve((zX, zY), label='Ground Return').opts(color='saddlebrown', tools=['hover'], height=550, width=400, line_width=2)\n", "ztVis = hv.Curve((zX, zT), label='RH100').opts(color='navy', tools=['hover'], height=550, width=400, line_width=2)\n", "rh25Vis = hv.Curve((zX, [rh25,rh25]),label='RH25').opts(color='lightblue',tools=['hover'], height=550, width=400, line_width=2)\n", "rh50Vis = hv.Curve((zX, [rh50,rh50]),label='RH50').opts(color='mediumblue',tools=['hover'], height=550, width=400, line_width=2)\n", "rh75Vis = hv.Curve((zX, [rh75,rh75]),label='RH75').opts(color='darkblue',tools=['hover'], height=550, width=400, line_width=2)" ] }, { "cell_type": "code", "execution_count": 44, "metadata": {}, "outputs": [ { "data": {}, "metadata": {}, "output_type": "display_data" }, { "data": { "application/vnd.holoviews_exec.v0+json": "", "text/html": [ "
\n", "
\n", "
\n", "" ], "text/plain": [ ":Overlay\n", " .Curve.Selected_Algorithm_left_parenthesis_a1_right_parenthesis :Curve [x] (y)\n", " .Curve.Ground_Return :Curve [x] (y)\n", " .Curve.RH100 :Curve [x] (y)\n", " .Curve.RH25 :Curve [x] (y)\n", " .Curve.RH50 :Curve [x] (y)\n", " .Curve.RH75 :Curve [x] (y)" ] }, "execution_count": 44, "metadata": { "application/vnd.holoviews_exec.v0+json": { "id": "beb3c2a5-0cb2-4119-a4ec-6919709febbe" } }, "output_type": "execute_result" } ], "source": [ "# Plot all of the metrics together\n", "l2aVis = rhVis * zVis * ztVis * rh25Vis * rh50Vis * rh75Vis\n", "l2aVis.opts(show_legend=True, legend_position='bottom_right', title='GEDI L2A Relative Height Metrics', ylabel='Elevation (m)',\n", " xlabel='Percent Energy Returned', xlim=(0, 100), ylim=(np.min(rhShot) + 1.5, np.max(rhShot) + 5), height=600,\n", " fontsize={'title':16, 'xlabel':16, 'ylabel': 16, 'legend': 14, 'xticks':12, 'yticks':12}, width=400)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Based on the figure above, it looks like the densest portion of the tree canopy for this shot is at around 45-55 meters. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 4.3 Combine RH Metrics and Waveforms" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Note that you will need to have downloaded the [waveform.csv](https://github.com/nasa/GEDI-Data-Resources/data/waveform.csv) file from the repo (if you have not cloned the repository) in order to execute the section below." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### In this section, we will import the corresponding waveform for our selected shot. " ] }, { "cell_type": "code", "execution_count": 45, "metadata": {}, "outputs": [], "source": [ "wvDF = pd.read_csv('data/waveform.csv')" ] }, { "cell_type": "code", "execution_count": 46, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Amplitude (DN)Elevation (m)
0227.00992111.252897
1226.57410111.103165
2226.66390110.953434
3227.26736110.803702
4228.11644110.653970
.........
1241230.74078-74.564040
1242230.04378-74.713771
1243229.24507-74.863503
1244228.71117-75.013235
1245228.53528-75.162966
\n", "

1246 rows × 2 columns

\n", "
" ], "text/plain": [ " Amplitude (DN) Elevation (m)\n", "0 227.00992 111.252897\n", "1 226.57410 111.103165\n", "2 226.66390 110.953434\n", "3 227.26736 110.803702\n", "4 228.11644 110.653970\n", "... ... ...\n", "1241 230.74078 -74.564040\n", "1242 230.04378 -74.713771\n", "1243 229.24507 -74.863503\n", "1244 228.71117 -75.013235\n", "1245 228.53528 -75.162966\n", "\n", "[1246 rows x 2 columns]" ] }, "execution_count": 46, "metadata": {}, "output_type": "execute_result" } ], "source": [ "wvDF" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Next, create a `holoviews` curve plot of the waveform. If you are interested in learning more about plotting L1B waveforms, be sure to check out the [GEDI L1B Tutorial](https://github.com/nasa/GEDI-Data-Resources/blob/main/python/tutorials/GEDI_L1B_V2_Tutorial.ipynb)." ] }, { "cell_type": "code", "execution_count": 47, "metadata": {}, "outputs": [ { "data": {}, "metadata": {}, "output_type": "display_data" }, { "data": { "application/vnd.holoviews_exec.v0+json": "", "text/html": [ "
\n", "
\n", "
\n", "" ], "text/plain": [ ":Curve [Amplitude (DN)] (Elevation (m))" ] }, "execution_count": 47, "metadata": { "application/vnd.holoviews_exec.v0+json": { "id": "710c5044-4e9f-4084-81db-834545a393ad" } }, "output_type": "execute_result" } ], "source": [ "# Create a holoviews interactive Curve plot with additional parameters defining the plot aesthetics \n", "visL1B = hv.Curve(wvDF).opts(color='darkgreen', tools=['hover'], height=600, width=400,\n", " xlim=(np.min(wvDF['Amplitude (DN)']) - 10, np.max(wvDF['Amplitude (DN)']) + 10), \n", " ylim=(np.min(wvDF['Elevation (m)']), np.max(wvDF['Elevation (m)'])),\n", " fontsize={'xticks':10, 'yticks':10,'xlabel':16, 'ylabel': 16, 'title':13}, line_width=2.5, title=f'{str(shot)}')\n", "visL1B" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Next, combine the L2A Relative Height (RH) metrics to the corresponding L1B waveform." ] }, { "cell_type": "code", "execution_count": 48, "metadata": {}, "outputs": [ { "data": {}, "metadata": {}, "output_type": "display_data" }, { "data": { "application/vnd.holoviews_exec.v0+json": "", "text/html": [ "
\n", "
\n", "
\n", "" ], "text/plain": [ ":Layout\n", " .Curve.I :Curve [Amplitude (DN)] (Elevation (m))\n", " .Overlay.I :Overlay\n", " .Curve.Selected_Algorithm_left_parenthesis_a1_right_parenthesis :Curve [x] (y)\n", " .Curve.Ground_Return :Curve [x] (y)\n", " .Curve.RH100 :Curve [x] (y)\n", " .Curve.RH25 :Curve [x] (y)\n", " .Curve.RH50 :Curve [x] (y)\n", " .Curve.RH75 :Curve [x] (y)" ] }, "execution_count": 48, "metadata": { "application/vnd.holoviews_exec.v0+json": { "id": "58119c4b-f58c-4200-8daf-427ef48ecb1c" } }, "output_type": "execute_result" } ], "source": [ "visL1B.opts(height=600, width=400, ylim=(np.min(rhShot), np.max(rhShot)+5), ylabel='Elevation (m)', xlabel='Amplitude (DN)') \\\n", "+ l2aVis.opts(height=600, width=400, ylim=(np.min(rhShot), np.max(rhShot)+5))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Above, notice the amplitude spike where the ground return is defined, and the close grouping of RH50 and RH75 over the dense portion of the upper canopy detected for this shot." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 4.4 Select Data from the Non-Default Algorithm\n", "#### Recall from above that the GEDI L2A product provides waveform processing results for multiple algorithm settings. Results for the selected algorithm selection are provided in the root directory of the data product for each beam. New with Version 2 is that the selected algorithm is selected as identifying the lowest non-noise mode (see Table 5 of the [GEDI L02 User Guide](https://lpdaac.usgs.gov/documents/988/GEDI02_User_Guide_V2.pdf)). Algorithm Setting Group 2 is the algorithm used for the data that we extracted in the sections above.\n", "> Elevation and height metrics outputs for all algorithm setting groups can be found in the `geolocation` subgroup of the L2A data product. \n", "For example: \n", " - elev_lowestreturn_a is the elevation of lowest return detected using algorithm setting group \n", " - rh_a are the relative height metrics at 1% intervals using algorithm (in cm). " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### In some cases, the selection of an alternative algorithm setting will provide a better result. See the [GEDI L02A User Guide](https://lpdaac.usgs.gov/documents/988/GEDI02_User_Guide_V2.pdf) for additional information." ] }, { "cell_type": "code", "execution_count": 49, "metadata": {}, "outputs": [], "source": [ "geolocationSDS = [b for b in beamSDS if '/geolocation/' in b] # Select all datasets within the geolocation subgroup" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Next, select the datasets for Algorithm Setting Group 3 (`a3`)." ] }, { "cell_type": "code", "execution_count": 50, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "['BEAM0110/geolocation/elev_highestreturn_a3',\n", " 'BEAM0110/geolocation/elev_lowestmode_a3',\n", " 'BEAM0110/geolocation/elev_lowestreturn_a3',\n", " 'BEAM0110/geolocation/elevs_allmodes_a3',\n", " 'BEAM0110/geolocation/energy_lowestmode_a3',\n", " 'BEAM0110/geolocation/lat_highestreturn_a3',\n", " 'BEAM0110/geolocation/lat_lowestmode_a3',\n", " 'BEAM0110/geolocation/lat_lowestreturn_a3',\n", " 'BEAM0110/geolocation/lats_allmodes_a3',\n", " 'BEAM0110/geolocation/lon_highestreturn_a3',\n", " 'BEAM0110/geolocation/lon_lowestmode_a3',\n", " 'BEAM0110/geolocation/lon_lowestreturn_a3',\n", " 'BEAM0110/geolocation/lons_allmodes_a3',\n", " 'BEAM0110/geolocation/num_detectedmodes_a3',\n", " 'BEAM0110/geolocation/quality_flag_a3',\n", " 'BEAM0110/geolocation/rh_a3',\n", " 'BEAM0110/geolocation/sensitivity_a3']" ] }, "execution_count": 50, "metadata": {}, "output_type": "execute_result" } ], "source": [ "a3SDS = [g for g in geolocationSDS if g.endswith('a3')] # Select algorithm 3 datasets\n", "a3SDS" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Repeat Section 4.2 for Algorithm Setting Group 3." ] }, { "cell_type": "code", "execution_count": 51, "metadata": {}, "outputs": [], "source": [ "# Bring in the desired SDS\n", "lata3 = gediL2A[[a for a in a3SDS if a.endswith('/lat_lowestmode_a3')][0]][index] # Latitude\n", "lona3 = gediL2A[[a for a in a3SDS if a.endswith('/lon_lowestmode_a3')][0]][index] # Latitude\n", "rhShot1 = gediL2A[[a for a in a3SDS if a.endswith('/rh_a3')][0]][index] / 100 # Relative height metrics\n", "zElevationa3 = gediL2A[[a for a in a3SDS if a.endswith('/elev_lowestmode_a3')][0]][index] # Elevation\n", "zTopa3 = gediL2A[[a for a in a3SDS if a.endswith('/elev_highestreturn_a3')][0]][index] # Elevation Highest Return\n", "rhShota3 = [z + zElevationa3 for z in rhShot1] # To convert canopy height to canopy elevation, add the elevation to each value" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Notice above that the non-default algorithm RH metrics are stored in cm, so we convert to meters by dividing by 100." ] }, { "cell_type": "code", "execution_count": 52, "metadata": {}, "outputs": [ { "data": {}, "metadata": {}, "output_type": "display_data" }, { "data": { "application/vnd.holoviews_exec.v0+json": "", "text/html": [ "
\n", "
\n", "
\n", "" ], "text/plain": [ ":Curve [x] (y)" ] }, "execution_count": 52, "metadata": { "application/vnd.holoviews_exec.v0+json": { "id": "b3ea791e-fdf8-4b83-a2e9-1ad7629a385a" } }, "output_type": "execute_result" } ], "source": [ "rhVisA3 = hv.Curve(rhShota3, label='Algorithm Setting Group 3 (a3)')\n", "rhVisA3 = rhVisA3.opts(color='green', tools=['hover'], height=500, width=400, title='GEDI L2A Relative Height Metrics (a3)', \n", " xlabel='Percent Energy Returned', ylabel='Elevation (m)', xlim=(0,100),ylim=(np.min(rhShot),np.max(rhShot)), \n", " fontsize={'title':14, 'xlabel':16, 'ylabel': 16, 'legend': 14, 'xticks':12, 'yticks':12}, line_width=3.5)\n", "rhVisA3" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Next, combine with the relative height metrics from the selected algorithm (Algorithm Setting Group 2):" ] }, { "cell_type": "code", "execution_count": 53, "metadata": {}, "outputs": [ { "data": {}, "metadata": {}, "output_type": "display_data" }, { "data": { "application/vnd.holoviews_exec.v0+json": "", "text/html": [ "
\n", "
\n", "
\n", "" ], "text/plain": [ ":Overlay\n", " .Curve.Selected_Algorithm_left_parenthesis_a1_right_parenthesis :Curve [x] (y)\n", " .Curve.Algorithm_Setting_Group_3_left_parenthesis_a3_right_parenthesis :Curve [x] (y)" ] }, "execution_count": 53, "metadata": { "application/vnd.holoviews_exec.v0+json": { "id": "16ecf5a5-bdc1-4abc-8bed-e264b28e878e" } }, "output_type": "execute_result" } ], "source": [ "(rhVis * rhVisA3).opts(show_legend=True, legend_position='bottom_right', title='GEDI L2A Relative Height Metrics by Algorithm',\n", " ylabel='Elevation (m)', xlabel='Percent Energy Returned', xlim=(0, 100), \n", " ylim=(np.min(rhShot) - 1.5, np.max(rhShot) + 1.5), height=600, width=600,\n", " fontsize={'title':16, 'xlabel':16, 'ylabel': 16, 'legend': 14, 'xticks':12, 'yticks':12})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Above we can see slight differences in the RH values between the two algorithms, particularly in the lower percentiles. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "# 5. Plot Transects " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Next, import a number of desired SDS layers for BEAM0110 (for the entire orbit) and create a `pandas` Dataframe to store the arrays." ] }, { "cell_type": "code", "execution_count": 54, "metadata": {}, "outputs": [], "source": [ "# Open all of the desired SDS\n", "dem = gediL2A[[g for g in beamSDS if g.endswith('/digital_elevation_model')][0]][()]\n", "srtm = gediL2A[[g for g in beamSDS if g.endswith('/digital_elevation_model_srtm')][0]][()]\n", "zElevation = gediL2A[[g for g in beamSDS if g.endswith('/elev_lowestmode')][0]][()]\n", "zHigh = gediL2A[[g for g in beamSDS if g.endswith('/elev_highestreturn')][0]][()]\n", "zLat = gediL2A[[g for g in beamSDS if g.endswith('/lat_lowestmode')][0]][()]\n", "zLon = gediL2A[[g for g in beamSDS if g.endswith('/lon_lowestmode')][0]][()]\n", "rh = gediL2A[[g for g in beamSDS if g.endswith('/rh')][0]][()]\n", "quality = gediL2A[[g for g in beamSDS if g.endswith('/quality_flag')][0]][()]\n", "degrade = gediL2A[[g for g in beamSDS if g.endswith('/degrade_flag')][0]][()]\n", "sensitivity = gediL2A[[g for g in beamSDS if g.endswith('/sensitivity')][0]][()]\n", "shotNums = gediL2A[f'{beamNames[0]}/shot_number'][()]\n", "selectedAlgorithm = gediL2A[[g for g in beamSDS if g.endswith('/selected_algorithm')][0]][()]\n", "\n", "# Create a shot index\n", "shotIndex = np.arange(shotNums.size)" ] }, { "cell_type": "code", "execution_count": 55, "metadata": {}, "outputs": [], "source": [ "canopyHeight = [r[100] for r in rh] # Grab RH100 (index 100 for each RH metrics)" ] }, { "cell_type": "code", "execution_count": 56, "metadata": {}, "outputs": [], "source": [ "# Take the DEM, GEDI-produced Elevation, and RH Metrics and add to a Pandas dataframe\n", "transectDF = pd.DataFrame({\n", " 'Shot Index': shotIndex,\n", " 'Shot Number': shotNums,\n", " 'Latitude': zLat,\n", " 'Longitude': zLon, \n", " 'Tandem-X DEM': dem,\n", " 'SRTM DEM': srtm,\n", " 'Elevation (m)': zElevation,\n", " 'Canopy Elevation (m)': zHigh, \n", " 'Canopy Height (rh100)': canopyHeight,\n", " 'Quality Flag': quality,\n", " 'Degrade Flag': degrade, \n", " 'Sensitivity': sensitivity,\n", " 'Selected Algorithm': selectedAlgorithm\n", "})" ] }, { "cell_type": "code", "execution_count": 57, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Shot IndexShot NumberLatitudeLongitudeTandem-X DEMSRTM DEMElevation (m)Canopy Elevation (m)Canopy Height (rh100)Quality FlagDegrade FlagSensitivitySelected Algorithm
002932060020041986926.923896-142.755692-999999.000000-999999.03693.5039063693.5039060.0005.0299311
112932060020041987026.924089-142.755501-999999.000000-999999.03693.5183113693.5183110.000-4.6206631
222932060020041987126.924283-142.755309-999999.000000-999999.03693.5327153693.5327150.0006.6137961
332932060020041987226.924477-142.755118-999999.000000-999999.03693.5471193693.5471190.000-0.7887311
442932060020041987326.924671-142.754927-999999.000000-999999.03693.5615233693.5615230.000-64.3275151
..........................................
1070091070092932060020052687851.797246-79.858241-36.021603-999999.04535.9052734535.9052730.00800.4461881
1070101070102932060020052687951.797246-79.857414-999999.000000-999999.04543.5742194543.5742190.00800.7098041
1070111070112932060020052688051.797245-79.856577-999999.000000-999999.04515.9028324515.9028320.0080-0.8169401
1070121070122932060020052688151.797245-79.855748-999999.000000-999999.04516.1020514516.1020510.0080-43.3016511
1070131070132932060020052688251.797245-79.854919-41.067085-999999.04515.5693364515.5693360.00800.4082801
\n", "

107014 rows × 13 columns

\n", "
" ], "text/plain": [ " Shot Index Shot Number Latitude Longitude Tandem-X DEM \\\n", "0 0 29320600200419869 26.923896 -142.755692 -999999.000000 \n", "1 1 29320600200419870 26.924089 -142.755501 -999999.000000 \n", "2 2 29320600200419871 26.924283 -142.755309 -999999.000000 \n", "3 3 29320600200419872 26.924477 -142.755118 -999999.000000 \n", "4 4 29320600200419873 26.924671 -142.754927 -999999.000000 \n", "... ... ... ... ... ... \n", "107009 107009 29320600200526878 51.797246 -79.858241 -36.021603 \n", "107010 107010 29320600200526879 51.797246 -79.857414 -999999.000000 \n", "107011 107011 29320600200526880 51.797245 -79.856577 -999999.000000 \n", "107012 107012 29320600200526881 51.797245 -79.855748 -999999.000000 \n", "107013 107013 29320600200526882 51.797245 -79.854919 -41.067085 \n", "\n", " SRTM DEM Elevation (m) Canopy Elevation (m) Canopy Height (rh100) \\\n", "0 -999999.0 3693.503906 3693.503906 0.0 \n", "1 -999999.0 3693.518311 3693.518311 0.0 \n", "2 -999999.0 3693.532715 3693.532715 0.0 \n", "3 -999999.0 3693.547119 3693.547119 0.0 \n", "4 -999999.0 3693.561523 3693.561523 0.0 \n", "... ... ... ... ... \n", "107009 -999999.0 4535.905273 4535.905273 0.0 \n", "107010 -999999.0 4543.574219 4543.574219 0.0 \n", "107011 -999999.0 4515.902832 4515.902832 0.0 \n", "107012 -999999.0 4516.102051 4516.102051 0.0 \n", "107013 -999999.0 4515.569336 4515.569336 0.0 \n", "\n", " Quality Flag Degrade Flag Sensitivity Selected Algorithm \n", "0 0 0 5.029931 1 \n", "1 0 0 -4.620663 1 \n", "2 0 0 6.613796 1 \n", "3 0 0 -0.788731 1 \n", "4 0 0 -64.327515 1 \n", "... ... ... ... ... \n", "107009 0 80 0.446188 1 \n", "107010 0 80 0.709804 1 \n", "107011 0 80 -0.816940 1 \n", "107012 0 80 -43.301651 1 \n", "107013 0 80 0.408280 1 \n", "\n", "[107014 rows x 13 columns]" ] }, "execution_count": 57, "metadata": {}, "output_type": "execute_result" } ], "source": [ "transectDF" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Notice the unusual values listed above--those shots are flagged as poor quality and will be removed in Section 5.1.\n", "#### Now that you have the desired SDS into a `pandas` dataframe, begin plotting the entire beam transect:" ] }, { "cell_type": "code", "execution_count": 58, "metadata": {}, "outputs": [ { "data": {}, "metadata": {}, "output_type": "display_data" }, { "data": { "application/vnd.holoviews_exec.v0+json": "", "text/html": [ "
\n", "
\n", "
\n", "" ], "text/plain": [ ":Scatter [x] (y)" ] }, "execution_count": 58, "metadata": { "application/vnd.holoviews_exec.v0+json": { "id": "3bce3113-ba53-4822-9293-d095e7916f48" } }, "output_type": "execute_result" } ], "source": [ "# Plot Canopy Height\n", "canopyVis = hv.Scatter((transectDF['Shot Index'], transectDF['Canopy Height (rh100)']))\n", "canopyVis.opts(color='darkgreen', height=500, width=900, title=f'GEDI L2A Full Transect {beamNames[0]}',\n", " fontsize={'title':16, 'xlabel':16, 'ylabel': 16}, size=0.1, xlabel='Shot Index', ylabel='Canopy Height (m)')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Congratulations! You have plotted your first GEDI sub-orbit beam transect. Notice above that things look a little messy--before we dive deeper into plotting full transects, let's quality filter the shots in the section below." ] }, { "cell_type": "code", "execution_count": 59, "metadata": {}, "outputs": [], "source": [ "del canopyVis, canopyHeight, degrade, dem, quality, sensitivity, shotIndex, shotNums, zElevation, zHigh, zLat, zLon" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 5.1 Quality Filtering\n", "#### Now that you have the desired layers imported as a dataframe for the entire beam transect, let's perform quality filtering.\n", "#### Below, remove any shots where the `quality_flag` is set to 0 by defining those shots as `nan`. \n", "#### The syntax of the line below can be read as: in the dataframe, find the rows \"where\" the quality flag is not equal (ne) to 0. If a row (shot) does not meet the condition, set all values equal to `nan` for that row." ] }, { "cell_type": "code", "execution_count": 60, "metadata": {}, "outputs": [], "source": [ "transectDF = transectDF.where(transectDF['Quality Flag'].ne(0)) # Set any poor quality returns to NaN" ] }, { "cell_type": "code", "execution_count": 61, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Shot IndexShot NumberLatitudeLongitudeTandem-X DEMSRTM DEMElevation (m)Canopy Elevation (m)Canopy Height (rh100)Quality FlagDegrade FlagSensitivitySelected Algorithm
0NaNNaNNaNNaNNaNNaNNaNNaNNaNNaNNaNNaNNaN
1NaNNaNNaNNaNNaNNaNNaNNaNNaNNaNNaNNaNNaN
2NaNNaNNaNNaNNaNNaNNaNNaNNaNNaNNaNNaNNaN
3NaNNaNNaNNaNNaNNaNNaNNaNNaNNaNNaNNaNNaN
4NaNNaNNaNNaNNaNNaNNaNNaNNaNNaNNaNNaNNaN
..........................................
107009NaNNaNNaNNaNNaNNaNNaNNaNNaNNaNNaNNaNNaN
107010NaNNaNNaNNaNNaNNaNNaNNaNNaNNaNNaNNaNNaN
107011NaNNaNNaNNaNNaNNaNNaNNaNNaNNaNNaNNaNNaN
107012NaNNaNNaNNaNNaNNaNNaNNaNNaNNaNNaNNaNNaN
107013NaNNaNNaNNaNNaNNaNNaNNaNNaNNaNNaNNaNNaN
\n", "

107014 rows × 13 columns

\n", "
" ], "text/plain": [ " Shot Index Shot Number Latitude Longitude Tandem-X DEM SRTM DEM \\\n", "0 NaN NaN NaN NaN NaN NaN \n", "1 NaN NaN NaN NaN NaN NaN \n", "2 NaN NaN NaN NaN NaN NaN \n", "3 NaN NaN NaN NaN NaN NaN \n", "4 NaN NaN NaN NaN NaN NaN \n", "... ... ... ... ... ... ... \n", "107009 NaN NaN NaN NaN NaN NaN \n", "107010 NaN NaN NaN NaN NaN NaN \n", "107011 NaN NaN NaN NaN NaN NaN \n", "107012 NaN NaN NaN NaN NaN NaN \n", "107013 NaN NaN NaN NaN NaN NaN \n", "\n", " Elevation (m) Canopy Elevation (m) Canopy Height (rh100) \\\n", "0 NaN NaN NaN \n", "1 NaN NaN NaN \n", "2 NaN NaN NaN \n", "3 NaN NaN NaN \n", "4 NaN NaN NaN \n", "... ... ... ... \n", "107009 NaN NaN NaN \n", "107010 NaN NaN NaN \n", "107011 NaN NaN NaN \n", "107012 NaN NaN NaN \n", "107013 NaN NaN NaN \n", "\n", " Quality Flag Degrade Flag Sensitivity Selected Algorithm \n", "0 NaN NaN NaN NaN \n", "1 NaN NaN NaN NaN \n", "2 NaN NaN NaN NaN \n", "3 NaN NaN NaN NaN \n", "4 NaN NaN NaN NaN \n", "... ... ... ... ... \n", "107009 NaN NaN NaN NaN \n", "107010 NaN NaN NaN NaN \n", "107011 NaN NaN NaN NaN \n", "107012 NaN NaN NaN NaN \n", "107013 NaN NaN NaN NaN \n", "\n", "[107014 rows x 13 columns]" ] }, "execution_count": 61, "metadata": {}, "output_type": "execute_result" } ], "source": [ "transectDF" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Below, quality filter even further by using the `degrade_flag` (Greater than zero if the shot occurs during a degrade period, zero otherwise) and the `Sensitivity` layer, using a threshold of 0.95." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### From the GEDI Level 2 User Guide:\n", "> Use the sensitivity metric available in L2A and L2B to select \"best\" data. The L2A and L2B quality_flag datasets use a conservative sensitivity threshold of 0.9 over land (0.5 over ocean), but under some conditions (e.g. dense forest) the user may benefit from selecting a higher threshold." ] }, { "cell_type": "code", "execution_count": 62, "metadata": {}, "outputs": [], "source": [ "transectDF = transectDF.where(transectDF['Degrade Flag'] < 1)\n", "transectDF = transectDF.where(transectDF['Sensitivity'] > 0.95)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Below, drop all of the shots that did not pass the quality filtering standards outlined above from the `transectDF`." ] }, { "cell_type": "code", "execution_count": 63, "metadata": {}, "outputs": [], "source": [ "transectDF = transectDF.dropna() # Drop all of the rows (shots) that did not pass the quality filtering above" ] }, { "cell_type": "code", "execution_count": 64, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Quality filtering complete, 15551 high quality shots remaining.\n" ] } ], "source": [ "print(f\"Quality filtering complete, {len(transectDF)} high quality shots remaining.\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 5.2 Plot Beam Transects\n", "#### Next, plot the full remaining transect of high quality values using `holoviews` Scatter(). Combine the Tandem-X derived elevation, SRTM elevation (New to Version 2!), the GEDI-derived elevation, and the Canopy Top Elevation in a combined holoviews plot." ] }, { "cell_type": "code", "execution_count": 65, "metadata": {}, "outputs": [], "source": [ "# Remove any shots where there is no Tandem-X Elevation values \n", "transectDF = transectDF.where(transectDF['Tandem-X DEM'] != -999999.0).dropna()\n", "\n", "# Plot Digital Elevation Model\n", "demVis = hv.Scatter((transectDF['Shot Index'], transectDF['Tandem-X DEM']), label='Tandem-X DEM')\n", "demVis = demVis.opts(color='black', height=500, width=900, fontsize={'xlabel':16, 'ylabel': 16}, size=1.5)" ] }, { "cell_type": "code", "execution_count": 66, "metadata": {}, "outputs": [], "source": [ "# Remove any shots where there is no SRTM Elevation values \n", "transectDF = transectDF.where(transectDF['SRTM DEM'] != -999999.0).dropna()\n", "\n", "# Plot SRTM Digital Elevation Model\n", "srtmVis = hv.Scatter((transectDF['Shot Index'], transectDF['SRTM DEM']), label='SRTM DEM')\n", "srtmVis = srtmVis.opts(color='darkblue', height=500, width=900, fontsize={'xlabel':16, 'ylabel': 16}, size=1.5)" ] }, { "cell_type": "code", "execution_count": 67, "metadata": {}, "outputs": [], "source": [ "# Plot GEDI-Retrieved Elevation\n", "zVis = hv.Scatter((transectDF['Shot Index'], transectDF['Elevation (m)']), label='GEDI-derived Elevation')\n", "zVis = zVis.opts(color='saddlebrown', height=500, width=900, fontsize={'xlabel':16, 'ylabel': 16}, size=1.5)" ] }, { "cell_type": "code", "execution_count": 68, "metadata": {}, "outputs": [], "source": [ "# Plot Canopy Top Elevation\n", "rhVis = hv.Scatter((transectDF['Shot Index'], transectDF['Canopy Elevation (m)']), label='Canopy Top Elevation')\n", "rhVis = rhVis.opts(color='darkgreen', height=500, width=900, fontsize={'xlabel':16, 'ylabel': 16}, size=1.5, \n", " tools=['hover'], xlabel='Shot Index', ylabel='Elevation (m)')" ] }, { "cell_type": "code", "execution_count": 69, "metadata": {}, "outputs": [ { "data": {}, "metadata": {}, "output_type": "display_data" }, { "data": { "application/vnd.holoviews_exec.v0+json": "", "text/html": [ "
\n", "
\n", "
\n", "" ], "text/plain": [ ":Overlay\n", " .Scatter.Tandem_hyphen_minus_X_DEM :Scatter [x] (y)\n", " .Scatter.SRTM_DEM :Scatter [x] (y)\n", " .Scatter.GEDI_hyphen_minus_derived_Elevation :Scatter [x] (y)\n", " .Scatter.Canopy_Top_Elevation :Scatter [x] (y)" ] }, "execution_count": 69, "metadata": { "application/vnd.holoviews_exec.v0+json": { "id": "1a906364-4858-46ee-8788-36d7618f58bc" } }, "output_type": "execute_result" } ], "source": [ "# Combine all three scatterplots\n", "(demVis * srtmVis * zVis * rhVis).opts(show_legend=True, legend_position='top_left',fontsize={'title':14, 'xlabel':16, \n", " 'ylabel': 16}, title=f'{beamNames[0]} Full Transect: {L2A.split(\".\")[0]}')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### The plot still looks a bit messy this far zoomed out--feel free to pan, zoom, and explore different areas of the plot. The shot plotted in section 4 was at index 45730. If you zoom into the high-quality shots between 4.000e+5 and 5.000e+5, you will find the portion of the transect intersecting Redwood National Park, seen below:\n", "![GEDI_L2A_V2_Tutorial_2.png](../../img/GEDI_L2A_V2_Tutorial_2.png \"BEAM0110 Transect Plot of Tandem-X Elevation, GEDI-derived Elevation, and Canopy Top Elevation over Redwood National Park, USA.\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 5.3 Subset Beam Transects" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Now, subset down to a smaller transect centered on the GEDI shot analyzed in the sections above." ] }, { "cell_type": "code", "execution_count": 70, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "45732\n" ] } ], "source": [ "print(index)" ] }, { "cell_type": "code", "execution_count": 71, "metadata": {}, "outputs": [], "source": [ "# Grab 50 points before and after the shot visualized above\n", "start = index - 50\n", "end = index + 50 " ] }, { "cell_type": "code", "execution_count": 72, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The transect begins at (41.27021331597735, -124.05764284824794) and ends at (41.30056975532455, -124.00258112950078).\n" ] } ], "source": [ "print(f\"The transect begins at ({transectDF['Latitude'][start]}, {transectDF['Longitude'][start]}) and ends at ({transectDF['Latitude'][end]}, {transectDF['Longitude'][end]}).\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Below, subset the transect using `.loc`." ] }, { "cell_type": "code", "execution_count": 73, "metadata": {}, "outputs": [], "source": [ "transectDF = transectDF.loc[start:end] # Subset the Dataframe to only the selected region of interest over Redwood NP" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 5.4 Plot RH Metrics Transects" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### In order to get an idea of the length of the beam transect that you are plotting, you can plot the x-axis as distance, which is calculated below." ] }, { "cell_type": "code", "execution_count": 74, "metadata": {}, "outputs": [], "source": [ "# Calculate along-track distance\n", "distance = np.arange(0.0, len(transectDF.index) * 60, 60) # GEDI Shots are spaced 60 m apart\n", "transectDF['Distance'] = distance # Add Distance as a new column in the dataframe" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### But how do you determine where the ground is? Next, create a similar plot but using the `rh` metrics from the L2A file. Notice that the `rh` metrics contain 100 values representing each percentile, so below we iterate through and append each value with the correct distance and elevation value to the list." ] }, { "cell_type": "code", "execution_count": 75, "metadata": {}, "outputs": [], "source": [ "# Create lists of canopy height, distance, and canopy elevation \n", "rhList, distList, ceList = [], [], []\n", "for s, i in enumerate(transectDF['Shot Index']):\n", " i = int(i)\n", " rhi = rh[i]\n", " for r in rhi:\n", " rhList.append(r)\n", " distList.append(distance[s])\n", " ceList.append(r + transectDF['Elevation (m)'][i]) # Convert to Canopy Top Elevation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Now plot canopy height vs. shot index as a scatter plot." ] }, { "cell_type": "code", "execution_count": 76, "metadata": {}, "outputs": [ { "data": {}, "metadata": {}, "output_type": "display_data" }, { "data": { "application/vnd.holoviews_exec.v0+json": "", "text/html": [ "
\n", "
\n", "
\n", "" ], "text/plain": [ ":Scatter [x] (y)" ] }, "execution_count": 76, "metadata": { "application/vnd.holoviews_exec.v0+json": { "id": "5c42fa93-dd93-4ea7-b04f-00e8c233660a" } }, "output_type": "execute_result" } ], "source": [ "title = 'Relative Height Metrics (0-100) across Redwood National Park: June 19, 2019'\n", "hv.Scatter((distList,rhList)).opts(color='darkgreen', alpha=0.1, height=500, width=900, size=7.5, \n", " xlabel='Distance Along Transect (m)', ylabel='Canopy Height (m)', title=title,\n", " fontsize={'title':16, 'xlabel':16, 'ylabel': 16, 'xticks':12, 'yticks':12})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Above, it looks like the transect is covering some tall trees! However, a quick google search confirms that Redwood trees can top out at over 100 m! The graph above is showing canopy *height*, but now add elevation to get a better idea of the overall three-dimensional structure of this transect across Redwood National Park." ] }, { "cell_type": "code", "execution_count": 77, "metadata": {}, "outputs": [], "source": [ "# Plot canopy elevation instead of canopy height\n", "rhVis = hv.Scatter((distList,ceList), label='RH 0-100').opts(color='darkgreen', alpha=0.1, height=500, width=900, size=5)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Next, add in the ground elevation." ] }, { "cell_type": "code", "execution_count": 78, "metadata": {}, "outputs": [ { "data": {}, "metadata": {}, "output_type": "display_data" }, { "data": { "application/vnd.holoviews_exec.v0+json": "", "text/html": [ "
\n", "
\n", "
\n", "" ], "text/plain": [ ":Overlay\n", " .Curve.Ground_Elevation :Curve [x] (y)\n", " .Scatter.RH_0_hyphen_minus_100 :Scatter [x] (y)" ] }, "execution_count": 78, "metadata": { "application/vnd.holoviews_exec.v0+json": { "id": "2b1dd603-80b4-4a07-8199-04847004db55" } }, "output_type": "execute_result" } ], "source": [ "# Plot GEDI-Retrieved Elevation\n", "demVis = hv.Curve((transectDF['Distance'], transectDF['Elevation (m)']), label='Ground Elevation')\n", "demVis = demVis.opts(color='black', height=500, width=900, line_width=3, xlabel='Distance Along Transect (m)',\n", " ylabel='Canopy Height (m)', fontsize={'title':14, 'xlabel':16, 'ylabel': 16, 'xticks':12, 'yticks':12},\n", " title='Relative Height Metrics and Elevation across Redwood National Park: June 19, 2019')\n", "\n", "# Plot Elevation and RH metrics together\n", "(demVis * rhVis).opts(legend_position='bottom_right')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Above, you can get an idea about the terrain over the region of interest, particularly the classic \"V\" representing the river valley that is bisected by the transect. In terms of vegetation structure, this plot does a good job of showing not only which portions of the canopy are taller, but also where they are denser (darker shades of green)." ] }, { "cell_type": "code", "execution_count": 79, "metadata": {}, "outputs": [], "source": [ "del ceList, distList, rhList, rhi" ] }, { "cell_type": "code", "execution_count": 80, "metadata": {}, "outputs": [], "source": [ "del beamSDS, distance, rh, transectDF" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### At this point you have visualized the elevation, canopy, and vertical structure of specific footprints over Redwood National Park, and for a transect cutting through the national park. In section 6 you will look at mapping all of the high-quality shots from all eight GEDI beams for a given region of interest in order to gain knowledge on the spatial distribution and characteristics of the canopy over Redwood National Park. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "# 6. Spatial Visualization\n", "#### Section 6 combines many of the techniques learned above including how to import GEDI datasets, perform quality filtering, spatial subsetting, and visualization. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 6.1 Import, Subset, and Quality Filter All Beams" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Below, re-open the GEDI L2A observation--but this time, loop through and import data for all 8 of the GEDI beams." ] }, { "cell_type": "code", "execution_count": 81, "metadata": {}, "outputs": [], "source": [ "beamNames = [g for g in gediL2A.keys() if g.startswith('BEAM')]" ] }, { "cell_type": "code", "execution_count": 82, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "['BEAM0000',\n", " 'BEAM0001',\n", " 'BEAM0010',\n", " 'BEAM0011',\n", " 'BEAM0101',\n", " 'BEAM0110',\n", " 'BEAM1000',\n", " 'BEAM1011']" ] }, "execution_count": 82, "metadata": {}, "output_type": "execute_result" } ], "source": [ "beamNames" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Loop through each of the desired datasets (SDS) for each beam, append to lists, and transform into a `pandas` DataFrame." ] }, { "cell_type": "code", "execution_count": 83, "metadata": {}, "outputs": [], "source": [ "# Set up lists to store data\n", "shotNum, dem, zElevation, zHigh, zLat, zLon, rh25, rh98, rh100 ,quality ,degrade, sensitivity ,beamI, selectedAlgorithm = ([] for i in range(14)) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Below, note that we are going to pull out the relative height metrics at the 25th, 98th, and 100th (canopy height) percentiles." ] }, { "cell_type": "code", "execution_count": 84, "metadata": {}, "outputs": [], "source": [ "# Loop through each beam and open the SDS needed\n", "for b in beamNames:\n", " [shotNum.append(h) for h in gediL2A[[g for g in gediSDS if g.endswith('/shot_number') and b in g][0]][()]]\n", " [dem.append(h) for h in gediL2A[[g for g in gediSDS if g.endswith('/digital_elevation_model') and b in g][0]][()]]\n", " [zElevation.append(h) for h in gediL2A[[g for g in gediSDS if g.endswith('/elev_lowestmode') and b in g][0]][()]] \n", " [zHigh.append(h) for h in gediL2A[[g for g in gediSDS if g.endswith('/elev_highestreturn') and b in g][0]][()]] \n", " [zLat.append(h) for h in gediL2A[[g for g in gediSDS if g.endswith('/lat_lowestmode') and b in g][0]][()]] \n", " [zLon.append(h) for h in gediL2A[[g for g in gediSDS if g.endswith('/lon_lowestmode') and b in g][0]][()]] \n", " [rh25.append(h[25]) for h in gediL2A[[g for g in gediSDS if g.endswith('/rh') and b in g][0]][()]] \n", " [rh98.append(h[98]) for h in gediL2A[[g for g in gediSDS if g.endswith('/rh') and b in g][0]][()]]\n", " [rh100.append(h[100]) for h in gediL2A[[g for g in gediSDS if g.endswith('/rh') and b in g][0]][()]] \n", " [quality.append(h) for h in gediL2A[[g for g in gediSDS if g.endswith('/quality_flag') and b in g][0]][()]] \n", " [degrade.append(h) for h in gediL2A[[g for g in gediSDS if g.endswith('/degrade_flag') and b in g][0]][()]] \n", " [sensitivity.append(h) for h in gediL2A[[g for g in gediSDS if g.endswith('/sensitivity') and b in g][0]][()]] \n", " [beamI.append(h) for h in [b] * len(gediL2A[[g for g in gediSDS if g.endswith('/shot_number') and b in g][0]][()])]\n", " [selectedAlgorithm.append(h) for h in gediL2A[[g for g in gediSDS if g.endswith('/selected_algorithm') and b in g][0]][()]]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Note that the cell above may take a few seconds to complete." ] }, { "cell_type": "code", "execution_count": 85, "metadata": {}, "outputs": [], "source": [ "# Convert lists to Pandas dataframe\n", "allDF = pd.DataFrame({\n", " 'Shot Number': shotNum,\n", " 'Beam': beamI,\n", " 'Latitude': zLat,\n", " 'Longitude': zLon,\n", " 'Tandem-X DEM': dem,\n", " 'Elevation (m)': zElevation,\n", " 'Canopy Elevation (m)': zHigh,\n", " 'Canopy Height (rh100)': rh100,\n", " 'RH 98': rh98,\n", " 'RH 25': rh25,\n", " 'Quality Flag': quality,\n", " 'Degrade Flag': degrade,\n", " 'Sensitivity': sensitivity,\n", " 'Selected Algorithm': selectedAlgorithm\n", "})" ] }, { "cell_type": "code", "execution_count": 86, "metadata": {}, "outputs": [], "source": [ "del beamI, degrade, dem, gediSDS, rh100, rh98, rh25, quality, sensitivity, zElevation, zHigh, zLat, zLon, shotNum, selectedAlgorithm" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 6.2 Spatial Subsetting\n", "#### Below, subset the pandas dataframe using a simple bounding box region of interest. If you are interested in spatially clipping GEDI shots to a GeoJSON region of interest, be sure to check out the GEDI-Subsetter python script available at: https://github.com/nasa/GEDI-Data-Resources/tree/main/python/scripts/GEDI_Subsetter." ] }, { "cell_type": "code", "execution_count": 87, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "790135" ] }, "execution_count": 87, "metadata": {}, "output_type": "execute_result" } ], "source": [ "len(allDF)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Over 3.5 million shots are contained in this single GEDI orbit! Below subset down to only the shots falling within this small bounding box encompassing Redwood National Park. `RedwoodNP` our `geopandas` geodataframe can be called for the \"envelope\" or smallest bounding box encompassing the entire region of interest. Here, use that as the bounding box for subsetting the GEDI shots." ] }, { "cell_type": "code", "execution_count": 88, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(-124.16015705494489,\n", " 41.080601363502545,\n", " -123.84950230520286,\n", " 41.83981133687605)" ] }, "execution_count": 88, "metadata": {}, "output_type": "execute_result" } ], "source": [ "redwoodNP.envelope[0].bounds" ] }, { "cell_type": "code", "execution_count": 89, "metadata": {}, "outputs": [], "source": [ "minLon, minLat, maxLon, maxLat = redwoodNP.envelope[0].bounds # Define the min/max lat/lon from the bounds of Redwood NP" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Filter by the bounding box, which is done similarly to filtering by quality in section 6.1 above." ] }, { "cell_type": "code", "execution_count": 90, "metadata": {}, "outputs": [], "source": [ "allDF = allDF.where(allDF['Latitude'] > minLat)\n", "allDF = allDF.where(allDF['Latitude'] < maxLat)\n", "allDF = allDF.where(allDF['Longitude'] > minLon)\n", "allDF = allDF.where(allDF['Longitude'] < maxLon)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### NOTE: It may take up to a few minutes to run the cell above." ] }, { "cell_type": "code", "execution_count": 91, "metadata": {}, "outputs": [], "source": [ "allDF = allDF.dropna() # Drop shots outside of the ROI" ] }, { "cell_type": "code", "execution_count": 92, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "4477" ] }, "execution_count": 92, "metadata": {}, "output_type": "execute_result" } ], "source": [ "len(allDF)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Notice you have drastically reduced the number of shots you are working with (which will greatly enhance your experience in plotting them below). But first, remove any poor quality shots that exist within the ROI." ] }, { "cell_type": "code", "execution_count": 93, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3079" ] }, "execution_count": 93, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Set any poor quality returns to NaN\n", "allDF = allDF.where(allDF['Quality Flag'].ne(0))\n", "allDF = allDF.where(allDF['Degrade Flag'] < 1) \n", "allDF = allDF.where(allDF['Sensitivity'] > 0.95)\n", "allDF = allDF.dropna()\n", "len(allDF)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Down to roughly 2000 shots, next create a `Shapely` Point out of each shot and insert it as the geometry column in the [soon to be geo]dataframe." ] }, { "cell_type": "code", "execution_count": 94, "metadata": {}, "outputs": [], "source": [ "# Take the lat/lon dataframe and convert each lat/lon to a shapely point\n", "allDF['geometry'] = allDF.apply(lambda row: Point(row.Longitude, row.Latitude), axis=1)" ] }, { "cell_type": "code", "execution_count": 95, "metadata": {}, "outputs": [], "source": [ "# Convert to geodataframe\n", "allDF = gp.GeoDataFrame(allDF)\n", "allDF = allDF.drop(columns=['Latitude','Longitude'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 6.3 Visualize All Beams: Canopy Height and Elevation" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Now, using the `pointVisual` function defined in section 3.2, plot the `geopandas` GeoDataFrame using `geoviews`." ] }, { "cell_type": "code", "execution_count": 96, "metadata": {}, "outputs": [ { "data": {}, "metadata": {}, "output_type": "display_data" }, { "data": { "application/vnd.holoviews_exec.v0+json": "", "text/html": [ "
\n", "
\n", "
\n", "" ], "text/plain": [ ":Overlay\n", " .WMTS.I :WMTS [Longitude,Latitude]\n", " .Points.I :Points [Longitude,Latitude] (Shot Number,Beam,Tandem-X DEM,Elevation (m),Canopy Elevation (m),Canopy Height (rh100),RH 98,RH 25,Quality Flag,Degrade Flag,Sensitivity,Selected Algorithm)\n", " .Polygons.I :Polygons [Longitude,Latitude]" ] }, "execution_count": 96, "metadata": { "application/vnd.holoviews_exec.v0+json": { "id": "f06d0a6c-6f33-4b99-88fe-29b22c233828" } }, "output_type": "execute_result" } ], "source": [ "allDF['Shot Number'] = allDF['Shot Number'].astype(str) # Convert shot number to string\n", "\n", "vdims = []\n", "for f in allDF:\n", " if f not in ['geometry']:\n", " vdims.append(f)\n", "\n", "visual = pointVisual(allDF, vdims = vdims)\n", "visual * gv.Polygons(redwoodNP['geometry']).opts(line_color='red', color=None)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Feel free to pan and zoom in to the GEDI shots in yellow. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Now let's plot the points in the geodataframe and add a colormap for Canopy Height (m) and Elevation (m). " ] }, { "cell_type": "code", "execution_count": 97, "metadata": {}, "outputs": [ { "data": {}, "metadata": {}, "output_type": "display_data" }, { "data": { "application/vnd.holoviews_exec.v0+json": "", "text/html": [ "
\n", "
\n", "
\n", "" ], "text/plain": [ ":Overlay\n", " .WMTS.I :WMTS [Longitude,Latitude]\n", " .Points.I :Points [Longitude,Latitude] (Shot Number,Beam,Tandem-X DEM,Elevation (m),Canopy Elevation (m),Canopy Height (rh100),RH 98,RH 25,Quality Flag,Degrade Flag,Sensitivity,Selected Algorithm)" ] }, "execution_count": 97, "metadata": { "application/vnd.holoviews_exec.v0+json": { "id": "0e28823d-78be-44db-8f49-c82a10976f48" } }, "output_type": "execute_result" } ], "source": [ "# Plot the basemap and geoviews Points, defining the color as the Canopy Height for each shot\n", "(gvts.EsriImagery * gv.Points(allDF, vdims=vdims).options(color='Canopy Height (rh100)',cmap='plasma', size=3, tools=['hover'],\n", " clim=(0,102), colorbar=True, clabel='Meters',\n", " title='GEDI Canopy Height over Redwood National Park: June 19, 2019',\n", " fontsize={'xticks': 10, 'yticks': 10, 'xlabel':16, 'clabel':12,\n", " 'cticks':10,'title':16,'ylabel':16})).options(height=500,\n", " width=900)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Above and in the screenshot below, notice the higher canopy heights (shades of yellow) over the Redwood stands of the national park vs. other types of forests (pink-blue) vs. the low-lying (and consequently flat) profiles over lakes and rivers (purple).\n", "![GEDI_L2A_V2_Tutorial_3.png](../../img/GEDI_L2A_V2_Tutorial_3.png \"GEDI Canopy Height over Redwood National Park: June 19, 2019.\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Next, take a look at the GEDI-derived elevation over the shots. Notice below that the colormap is changed to 'terrain'. " ] }, { "cell_type": "code", "execution_count": 98, "metadata": {}, "outputs": [ { "data": {}, "metadata": {}, "output_type": "display_data" }, { "data": { "application/vnd.holoviews_exec.v0+json": "", "text/html": [ "
\n", "
\n", "
\n", "" ], "text/plain": [ ":Overlay\n", " .WMTS.I :WMTS [Longitude,Latitude]\n", " .Points.I :Points [Longitude,Latitude] (Shot Number,Beam,Tandem-X DEM,Elevation (m),Canopy Elevation (m),Canopy Height (rh100),RH 98,RH 25,Quality Flag,Degrade Flag,Sensitivity,Selected Algorithm)" ] }, "execution_count": 98, "metadata": { "application/vnd.holoviews_exec.v0+json": { "id": "0966dde7-7874-41bc-bc9d-625e4caadecc" } }, "output_type": "execute_result" } ], "source": [ "(gvts.EsriImagery * gv.Points(allDF, vdims=vdims).options(color='Elevation (m)',cmap='terrain', size=3, tools=['hover'],\n", " clim=(min(allDF['Elevation (m)']), max(allDF['Elevation (m)'])),\n", " colorbar=True, clabel='Meters',\n", " title='GEDI Elevation over Redwood National Park: June 19, 2019',\n", " fontsize={'xticks': 10, 'yticks': 10, 'xlabel':16, 'clabel':12,\n", " 'cticks':10,'title':16,'ylabel':16})).options(height=500,\n", " width=900)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Success! You have now learned how to start working with Version 2 GEDI L2A files in Python as well as some interesting strategies for visualizing those data in order to better understand your specific region of interest. Using this Jupyter Notebook as a workflow, you should now be able to switch to GEDI files over your specific region of interest and re-run the notebook. Good Luck!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "# 7. Export Subsets as GeoJSON Files\n", "#### In this section, export the GeoDataFrame as a `.geojson` file that can be easily opened in your favorite remote sensing and/or GIS software and will include an attribute table with all of the shots/values for each of the SDS layers in the dataframe." ] }, { "cell_type": "code", "execution_count": 99, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'data/GEDI02_A_2019170155833_O02932_02_T02267_02_003_01_V002.h5'" ] }, "execution_count": 99, "metadata": {}, "output_type": "execute_result" } ], "source": [ "gediL2A.filename # L2A Filename" ] }, { "cell_type": "code", "execution_count": 100, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'data/GEDI02_A_2019170155833_O02932_02_T02267_02_003_01_V002.json'" ] }, "execution_count": 100, "metadata": {}, "output_type": "execute_result" } ], "source": [ "outName = gediL2A.filename.replace('.h5', '.json') # Create an output file name using the input file name\n", "outName" ] }, { "cell_type": "code", "execution_count": 101, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "c:\\Users\\mjami\\AppData\\Local\\miniforge3\\envs\\lpdaac_vitals\\Lib\\site-packages\\pyogrio\\geopandas.py:662: UserWarning: 'crs' was not provided. The output dataset will not have projection information defined and may not be usable in other systems.\n", " write(\n" ] } ], "source": [ "allDF.to_file(outName, driver='GeoJSON') # Export to GeoJSON" ] }, { "cell_type": "code", "execution_count": 102, "metadata": {}, "outputs": [], "source": [ "del allDF " ] }, { "cell_type": "markdown", "metadata": { "jp-MarkdownHeadingCollapsed": true, "tags": [] }, "source": [ "\n", "## Contact Info: \n", "\n", "Email: LPDAAC@usgs.gov \n", "Voice: +1-866-573-3222 \n", "Organization: Land Processes Distributed Active Archive Center (LP DAAC)¹ \n", "Website: \n", "Date last modified: 02-20-2024 \n", "\n", "¹Work performed under USGS contract G15PD00467 for NASA contract NNG14HH33I. " ] } ], "metadata": { "kernelspec": { "display_name": "lpdaac_vitals", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.10" } }, "nbformat": 4, "nbformat_minor": 4 }