# Wyvern Knowledge Centre — Full Text > The Wyvern Knowledge Centre teaches how to work with Wyvern hyperspectral satellite Earth observation imagery. It contains background articles on GIS and hyperspectral concepts, hands-on tutorials for QGIS/EnMAP-Box, Python, ArcGIS Pro, and ENVI, real-world use cases, and an interactive hyperspectral index library. All tutorials use freely available imagery from Wyvern's Open Data Program. This file concatenates the full text of every page in the Knowledge Centre. Source: https://knowledge.wyvern.space --- [Source: https://knowledge.wyvern.space/docs/tutorials/arcgis_pro/] # ArcGIS Pro [ArcGIS Pro](https://pro.arcgis.com/en/pro-app/latest/get-started/get-started.htm) is a full-featured GIS software that allows you to explore, visualize, and analyze data. ### Tutorials ### More Resources ArcGIS Pro has extensive [documentation and tutorials](https://www.esri.com/en-us/arcgis/products/arcgis-pro/resources). --- [Source: https://knowledge.wyvern.space/docs/tutorials/envi/] # ENVI [ENVI](https://www.nv5geospatialsoftware.com/Products/ENVI) is a fully featured GIS software package with a focus on processing & analyzing hyperspectral imagery. ## Tutorials - [🛰️ Using Dragonette Imagery Data In ENVI Software](dragonette_data_in_envi): Data ingest, display, atmospheric correction, and spectral analysis of Dragonette imagery in ENVI ## More Resources NV5 provides [excellent documentation](https://www.nv5geospatialsoftware.com/docs/using_envi_Home.html) for their ENVI software package. There is also fantastic documentation around [different vegetative indices](https://www.nv5geospatialsoftware.com/docs/NarrowbandGreenness.html) & [processing algorithms](https://www.nv5geospatialsoftware.com/docs/AnalyzeData.html). --- [Source: https://knowledge.wyvern.space/docs/tutorials/python/] # Python Python is a fantastic tool for analyzing geospatial data and building applications. It has an incredibly rich and active community for all topics in computing. ## Tutorials ## Relevant Python Packages Here we'll focus on some popular and useful packages for processing geospatial data. There is also a huge wealth of other packages out there waiting to be used! Don't be afraid to search for packages tailored to your specific use case. ### Rasterio & GDAL :::tip GDAL can be difficult to install on Windows systems! One easy way around this problem is using [Anaconda](https://www.anaconda.com/) instead of the regular Python package system. ::: [Rasterio](https://rasterio.readthedocs.io/en/latest/index.html) is a Python library for reading, writing, and manipulating geospatial raster data. It is built on top of the popular GDAL library, which is a powerful tool for working with geospatial data of all kinds. Rasterio makes it easy to read and write raster data in various formats, and provides a range of functions and methods for manipulating and analyzing raster data. You can use Rasterio to perform operations like cropping, resampling, and reprojecting raster data, as well as to extract metadata and other information from raster datasets. Rasterio is a valuable tool for anyone working with geospatial raster data, whether you are a GIS specialist, data scientist, or just want to perform some simple raster analysis tasks. ### NumPy [Numpy](https://numpy.org/doc/stable/index.html) is a powerful Python library for working with large, multi-dimensional arrays and matrices of numerical data (exactly what geospatial rasters are!). It is a key component of many scientific computing projects and is widely used in machine learning and data analysis. Numpy provides a range of functions and methods for creating and manipulating arrays, as well as performing mathematical operations on them. It is known for its efficiency and performance, making it a popular choice for working with large datasets. :::info Numpy is able to compute complex array operations like index calculations for NDVI, NDWI, VREI2, etc! Check out the article on [Index Calculation](index_calculation) and [Plant Health Indices](/docs/documentation/datasets_and_use_cases/plant_health_indices) for more details on calculating indices. ::: ### GeoPandas [Geopandas](https://geopandas.org/en/stable/) is a Python library that allows you to work with geospatial data, such as maps and geographical information. It is built on top of the popular data manipulation library pandas, and makes it easy to work with geospatial data in a pandas-like way. Geopandas allows you to perform spatial operations, such as merging and joining spatial data (vector data especially), as well as visualize and plot spatial data on maps. Geopandas is built on top of other popular libraries like pandas and numpy, and is an essential part of many data analysis projects. ### Xarray [Xarray](https://docs.xarray.dev/en/stable/) is a Python library for working with multi-dimensional arrays and labeled data. It is particularly useful for working with large and complex datasets, such as those found in the Earth sciences, and offers a range of tools for manipulating and analyzing such data. Xarray introduces labeled dimensions and coordinates to NumPy arrays, which allows you to work with data that has explicit meaning beyond just its numerical values. This makes it easier to understand and work with your data, and also enables you to perform operations that preserve the metadata and structure of your data. Xarray is built on top of NumPy and other scientific computing libraries, and is a valuable tool for anyone working with multi-dimensional data. ### Dask [Dask](https://www.dask.org/) is a flexible parallel computing library for Python that allows you to scale your computations beyond the limits of a single machine. It is particularly useful for working with large and complex datasets that don't fit in memory, and can be used as a drop-in replacement for NumPy, pandas, or other popular scientific computing libraries. Dask provides a range of tools for parallelizing your code and distributing your computations across multiple cores or even distributed across a cluster of machines. It is a powerful tool for data scientists and analysts working with large datasets, and is often used in combination with other libraries like NumPy and pandas. ### Jupyter Notebooks [Jupyter notebooks](https://jupyter.org/) are useful for hyperspectral analysis with Python because they provide a flexible platform for working with geospatial data and libraries. For example, notebooks can be used to load, manipulate, and visualize geospatial data using libraries such as GeoPandas and Folium. This allows users to quickly explore and analyze their data in a more interactive and intuitive way. Additionally, notebooks support the ability to mix code, text, and visualizations in a single document, making it easy to communicate and share geospatial analysis and results. Overall, Python notebooks are a powerful tool for working with geospatial data and performing geospatial analysis. ### Spectral Python (SPy) [Spectral Python](https://www.spectralpython.net/) is a Python library for interacting with hyperspectral imagery. It is able to read in a number of file formats, including ENVI and AVIRIS. SPy has functions for a number of different algorithms, including dimensionality reduction, and classification using supervised & unsupervised methods. SPy is also capable of reading in spectral libraries for materials, including building ENVI spectral libraries. ### More Resources [python-geospatial (list compiled by Qiusheng Wu)](https://github.com/giswqs/python-geospatial) --- [Source: https://knowledge.wyvern.space/docs/tutorials/qgis/] # QGIS [QGIS](https://www.qgis.org/) is a free, fully featured, and open-source geospatial software that allows you to perform geospatial analysis on a wide variety of data and topics. It is compatible with Windows, Linux, and Mac. [EnMAP-Box](https://enmap-box.readthedocs.io/en/latest/general/about.html) is a free and open source plugin for QGIS, developed at Humboldt-Universität zu Berlin and Universität Greifswald. It is focused on processing hyperspectral/imaging spectroscopy data, especially data produced by the [EnMAP](https://www.enmap.org/) satellite. Many of the tools and functions in EnMAP-Box are also compatible with Wyvern's imagery! ### Tutorials ### More Resources EnMAP has incredibly fantastic [documentation](https://enmap-box.readthedocs.io/en/latest/usr_section/usr_gettingstarted.html) and [education materials](https://www.enmap.org/events_education/hyperedu). --- [Source: https://knowledge.wyvern.space/docs/documentation/datasets_and_use_cases/plant_health_indices] # Plant health indices A broadband index like NDVI tells you how much green vegetation is in a pixel. A narrowband index can tell you what state it's in: how much chlorophyll the leaves hold, whether stress pigments are building up, and where the red edge sits. This page explains why that difference exists, which plant properties you can measure, and which indices Dragonette supports. ## Background Multispectral satellites such as Landsat and Sentinel-2 are the usual source for vegetation indices like [NDVI](https://en.wikipedia.org/wiki/Normalized_difference_vegetation_index), [SAVI](https://en.wikipedia.org/wiki/Soil-adjusted_vegetation_index), and [EVI](https://en.wikipedia.org/wiki/Enhanced_vegetation_index). These measure the presence and general vigour of vegetation. What they can't do reliably is separate the individual things that make up plant health, such as chlorophyll content, carotenoid content, or leaf area, because the absorption features of those pigments are narrower than the bands and overlap with one another. NDVI still does its own job well. In the images below it separates vegetation cleanly from everything else, picks out areas of denser growth, and even shows crop rows within fields. Sentinel-2 RGB Sentinel-2 NDVI Sentinel-2 true colour image and the NDVI calculated from it Figure 1 shows why. A near-infrared band separates green vegetation from soil easily, because vegetation reflects far more near-infrared light. But the bands are wide: Sentinel-2's Band 8, the near-infrared band marked 4 in the figure, is about 100 nm across. Figure 1: Typical multispectral bands (grey) against soil, vegetation, and water spectra. Each band averages over a wide stretch of the spectrum. `[2]` Figure 2 shows how much detail a band that wide passes over. It marks the wavelengths researchers link to specific plant properties, and in the visible and red edge many of them are only a few tens of nanometres apart: pigments at 450, 515, 550, 570, and 650 nm, chlorophyll absorption at 687 nm, and stress at 705 and 720 nm. A 100 nm band averages several of them into one number. Narrow bands measure them separately. Figure 2: Reflectance of four crops at different growth stages, annotated with the wavelengths linked to specific plant properties. Dragonette covers the region from 445 to 870 nm. `[3]` ## What you can measure Figure 3: Absorbance of chlorophyll b, carotene, and anthocyanin. Each peaks in a different place, but the curves overlap, so separating them takes bands narrow enough to sit on each peak. `[1]` These are the properties of foliage that reflectance can measure: - [Chlorophyll](https://en.wikipedia.org/wiki/Chlorophyll): the green pigments a plant uses to absorb light for photosynthesis. Chlorophyll content is the most direct measure of how actively a plant is growing. - [Carotenoids](https://en.wikipedia.org/wiki/Carotenoid): yellow, orange, and red pigments with several roles, including protecting leaves from excess light. A high ratio of carotenoids to chlorophyll is a sign of stress. - [Anthocyanins](https://en.wikipedia.org/wiki/Anthocyanin): red, purple, and blue pigments that are mostly absent through the growing season and build up towards its end. Along with carotenoids, they give leaves their autumn colours. - [Leaf area index (LAI)](https://en.wikipedia.org/wiki/Leaf_area_index): the area of green leaf per unit of ground area. A dimensionless measure of how dense the canopy is. - [Red edge position](https://en.wikipedia.org/wiki/Red_edge): the wavelength of the steep rise in reflectance between red and near-infrared, roughly 680 to 750 nm. It shifts with chlorophyll content, which makes it useful for crop identification, growth stage, disturbance detection, and canopy stress. - [Senescence](https://en.wikipedia.org/wiki/Plant_senescence#:~:text=Plant%20senescence%20is%20the%20process,leaf%20color%20in%20deciduous%20trees.): the aging of a plant or leaf, usually marked by falling chlorophyll and rising anthocyanins. Each of these has indices designed to measure it. The table lists popular ones and which satellites can compute them. :::note This is a small selection. The [Index Library](/hyperspectral_library) has the full set, with formulas and per-instrument band mappings. ::: | Index | Applications | Sentinel‑2 | Landsat 7/8 | SuperDove | Wyvern Dragonette | | ---- | ---- | :----: | :----: | :----: | :----: | | **NDVI** | General plant health & presence of vegetation | | | | * | | **EVI** | General plant health & presence of vegetation, counteracts NDVI saturation problem | | | | * | | **SAVI** | General plant health & presence of vegetation, suppresses effect of soil pixels | | | | * | | **RENDVI** | Sensitive to small changes in senescence & foliage content | | | | | | **VREI2** | Sensitive to chlorophyll concentration, leaf area index, and water content | | | | | | **REPI** | Crop identification, disturbance detection, canopy stress | | | | | | **ARI1/2** | Sensitive to weakening vegetation due to high concentration of anthocyanins | | | | | | **CRI1/2** | Sensitive to stressed vegetation due to high concentration of carotenoids | | | | | | **SIPI/1** | Sensitive to ratio of carotenoids to chlorophyll | | | | | ## Reading the table The indices above fall into three groups, and knowing which group you are in tells you what the number can and cannot say. **Broadband greenness** (NDVI, EVI, SAVI) measures how much green vegetation is present. Every sensor in the table can compute it, including Dragonette, which is what the asterisk marks: these indices were designed around wide bands, so on a hyperspectral sensor you either average several narrow bands together to imitate the original broadband response, or you substitute a narrowband index that measures the same thing better. Both approaches work, but neither one takes full advantage of hyperspectral data. **Red edge indices** (RENDVI, VREI2, REPI) measure where the reflectance climb between red and near-infrared sits, and how steep it is. That position shifts with chlorophyll concentration and leaf area, and it shifts early, before a canopy looks any different to the eye or to a broadband sensor. Dragonette carries seven bands between 690 and 765 nm, so the red edge is sampled repeatedly across its width rather than straddled by one wide band. **Pigment indices** (ARI1/2, CRI1/2, SIPI/1) target individual compounds: anthocyanins, carotenoids, and the ratio of carotenoids to chlorophyll. Their absorption features are narrow and close together in the 500 to 800 nm range, so they need bands placed precisely on the feature. A wide band that spans the feature and the flat spectrum either side averages the signal away. The pattern down the Wyvern column follows from that. It is not that broadband satellites are poor instruments; it is that an index can only be computed if the sensor has bands at the wavelengths the formula uses, and narrow bands mean the formula can be applied where the feature actually is instead of near it. Two caveats worth carrying into any of these. First, an index is a ratio, not a measurement: it responds to the feature you are after, and also to soil background, canopy geometry, illumination, and residual atmosphere. Compare values within a scene before comparing across dates. Second, all of these want L2A surface reflectance. On L1B radiance the solar spectrum and atmospheric absorption are larger than the differences the index is measuring, so the result is not meaningful. See [processing levels](../getting_started/processing_levels) for the difference. ## Try it yourself Formulas, band mappings for each Dragonette configuration, and the source papers for every index above are in the [Index Library](/hyperspectral_library). To compute one on a real scene, the [index calculation tutorial](/docs/tutorials/python/index_calculation) walks through RENDVI on a Dragonette L2A scene in Python, selecting bands by wavelength rather than by position. The [water quality indices tutorial](/docs/tutorials/python/water_quality_indices) applies the same approach to a different problem. ## Citations `[1]` Zhou, Ning & Dong, Weiming & Mei, Xing. (2006). Realistic Simulation of Seasonal Variant Maples. 295 - 301. 10.1109/PMA.2006.29. `[2]` Pham, Binh. (2018). Satellite remote sensing of the variability of the continental hydrology cycle in the lower Mekong basin over the last two decades. 10.13140/RG.2.2.27136.07686. `[3]` [Hyperspectral Signatures by USGS Western Geographic Science Center](https://www.usgs.gov/media/images/hyperspectral-signatures-0) --- [Source: https://knowledge.wyvern.space/docs/documentation/getting_started/geospatial_basics] # Geospatial basics Wyvern scenes are rasters. Your field boundaries, sample points, and masks are vectors. Almost every analysis you'll do combines the two, so it's worth being clear on how they differ before moving on to hyperspectral imagery. ## Vector data Vector data represents features as points, lines, and polygons, each stored as a list of coordinates. Roads, buildings, field boundaries, and sampling sites are all vectors. Because a vector feature is geometry rather than a grid of cells, it stays sharp at any zoom level, and it's small and quick to edit, reproject, and query. Each feature can also carry attributes, such as a field's crop type or a site's sample ID, which is what you colour or filter a map by. Common formats are GeoJSON, GeoPackage, and shapefile. Vector data: grazing area polygons in Alberta, labelled from each polygon's name attribute. ## Raster data Raster data is a grid of cells, or pixels, each holding a value for one location on the ground. Elevation models, land cover maps, temperature grids, and every satellite image are rasters. A raster can have several layers, called bands, stacked over the same grid. A Dragonette scene has 23 or 31 of them, one per wavelength. Rasters are larger than vectors and slower to process, and their detail is fixed by the pixel size. Zooming in on a 5.3 m pixel doesn't reveal anything smaller than 5.3 m. It just makes the square bigger. Common formats are GeoTIFF, including the Cloud-Optimized GeoTIFF (COG) that Wyvern delivers, NetCDF, and HDF5. See [hyperspectral image formats](hyperspectral_image_formats) for more. Raster data: a true colour image of farmland in Germany. Every field boundary you can see is a pattern in the pixels, not a stored shape. ## Using them together Most analyses use vectors to ask questions of rasters. You clip a scene to a field boundary, average an index inside each polygon, or pull a spectrum at each sample point. You can also go the other way and turn a raster result into vectors, for example by outlining every area where an index crosses a threshold. The [water quality tutorial](/docs/tutorials/python/water_quality_indices) ends this way, exporting areas of high chlorophyll as a GeoJSON file you can load into any GIS. Both kinds of data need to share a coordinate reference system before you combine them. Wyvern L2A scenes are delivered in the UTM zone that contains the scene, while most vector data you download will be in latitude and longitude (EPSG:4326). QGIS reprojects on the fly for display, but for analysis in code, reproject one to match the other first. Vector crop-type polygons overlaid on a raster image. This pairing is how you get labelled training data for a classifier, or per-field statistics from an index. ## More reading [Intro to GIS and Spatial Analysis by Manuel Gimond](https://mgimond.github.io/Spatial/index.html) [Colorado State GIS Basics](https://gis.colostate.edu/resources/learning-resources/) --- [Source: https://knowledge.wyvern.space/docs/documentation/more/material_reflectance] # Material reflectance libraries Every material reflects light differently depending on its chemical composition and structure. A spectral library is a collection of those reflectance spectra, measured on known samples, and it's what you compare image spectra against to identify materials. See [hyperspectral signatures](../getting_started/hyperspectral_signatures) for how that comparison works, and remember to resample library spectra to Dragonette's bands before you compare. ## USGS Spectral Library The United States Geological Survey's [Spectral Library Version 7](https://crustal.usgs.gov/speclab/QueryAll07a.php) (splib07) holds spectra for close to 2,500 samples, from the ultraviolet to the shortwave infrared. Minerals are its core, and it also covers soils, vegetation, liquids, organic compounds including hydrocarbons, and artificial materials. You can search and download spectra on the site. The USGS Spectral Library Version 7 search page. ## ECOSTRESS Spectral Library [The ECOSTRESS Spectral Library](https://speclib.jpl.nasa.gov/) from NASA JPL contains over 3,400 spectra of natural and human-made materials. It incorporates three earlier libraries, from Johns Hopkins University, JPL, and the USGS. ## Urban surfaces spectral library The University of Reading's London Urban Meteorological Observatory has [published a spectral library](https://micromet.reading.ac.uk/spectral-library/) focused on urban materials. ## OpenSpecLib [OpenSpecLib](https://github.com/null-jones/openspeclib) is an open source project that packages the USGS, ECOSTRESS, and EcoSIS libraries as Parquet files, so you can load them straight into Python. The [rare earth elements tutorial](/docs/tutorials/python/detecting_rare_earth_elements) uses it to pull a USGS reference spectrum and resample it to Dragonette's bands. --- [Source: https://knowledge.wyvern.space/docs/documentation/] # Knowledge Centre Welcome to Wyvern's Knowledge Centre 🌈. This section explains the ideas you need before you open a hyperspectral scene: how geospatial data is structured, what hyperspectral imaging measures, how our data is processed and delivered, and what you can do with it. When you're ready to work with real data, the [tutorials](/docs/tutorials/) pick up from here. :::info Questions, corrections, or something you'd like us to write about? Post on [GitHub Discussions](https://github.com/Nrevyw/wyvern-public-resources/discussions) or find us on [LinkedIn](https://www.linkedin.com/company/wyvern-space). ::: --- [Source: https://knowledge.wyvern.space/docs/tutorials/arcgis_pro/loading_data] # Loading Data in ArcGIS Pro Final output of this tutorial ## Before you start... Before starting this tutorial, please ensure you have done the following: - ArcGIS Pro has been downloaded installed on your machine. Visit [pro.arcgis.com](https://pro.arcgis.com/en/pro-app/latest/get-started/download-arcgis-pro.htm) to download ArcGIS Pro. - You have a Wyvern image you can load into ArcGIS Pro. If you haven't downloaded one already, we'll be [using this image](https://opendata.wyvern.space/#/year/2024/wyvern_dragonette-001_20241003T001114_b9b1a0b8/wyvern_dragonette-001_20241003T001114_b9b1a0b8.json) for the next steps. In this tutorial, there will be references to the [Ribbon](https://www.esri.com/arcgis-blog/products/arcgis-pro/uncategorized/arcgis-pro-ribbons-toolbars-and-ui-hacks/). This is the primary interface for working with ArcGIS Pro. ArcGIS Ribbon ## Starting a new ArcGIS Pro Project Before opening the Wyvern image, an ArcGIS Project should be created: - Open ArcGIS Pro. - Click `New Project` -> `Map`. - Name your project and chose a location for saved files on your computer. ArcGIS Pro will now be open with the basemaps `World Topographic Map` and `World Hillshade`. Empty ArcGIS Pro Project ## Loading Data into ArcGIS Pro To load an image into ArcGIS, you have a number of different options: - Using the `Add Data` tool. - `Ribbon` -> `Map` -> `Add Data` -> `Browse` -> Navigate to and select the Wyvern GeoTIFF file -> `Click OK`. - From your file explorer. - Drag-and-drop the GeoTIFF file into the ArcGIS Pro main window. - Depending on your ArcGIS settings, you may be prompted to `Calculate Statistics`. This option calculates image statistics which will allow ArcGIS to more quickly read and display the image. It will take some time to calculate once but will be faster in the long run! Congrats! You have successfully loaded an image into ArcGIS Pro. ArcGIS will automatically read wavelength information about the Dragonette image. This means it will automatically load the correct bands for an natural color RGB! Wyvern GeoTIFF rendered in ArcGIS Pro Now, to see more detail, let's zoom into an area we're interested in using the scroll wheel or zoom tools. Zoom of Wyvern hyperspectral data natural color ### Band Combinations Since ArcGIS automatically reads Wyvern's band information, it can also quickly change between natural color and color infrared band combinations! Below is a list of bands that ArcGIS will use for visualizations. | Name | Description | Red Band | Green Band | Blue Band | | ---- | ----------- | -------- | ---------- | --------- | | Natural Color (RGB) | Normal visualization | 660nm | 549nm | 503nm | | Color Infrared | Vegetation analysis | 799nm | 660nm | 549nm | To change the band combinations to color infrared, first select the Wyvern image in the `Contents Pane` to the left of the map. Then click `Ribbon` -> `Raster Layer` -> `Band Combination` -> `Color Infrared`. Color Infrared Image ### Adjusting Stretch Every image is slightly different in its lighting conditions and distribution of values within each band. ArcGIS has default stretch settings that work alright for most images, but definitely benefit from careful adjustment. To change the stretch, select the image layer in the `Contents Pane` to the left of the map. Then `right-click the image layer` -> `Symbology`. The symbology pane will open to the right of the map. Opening the Symbology Pane In this case, the `Stretch type` is set to `Percent Clip`. To easily adjust the stretch, change `Stretch type` to `Standard Deviation`. This will default the `Number of standard deviations` to `2`. This value can be changed to increase or decrease the difference between bright and dark pixels. Stretch standard deviation 2 In this image, the bright clouds are affecting the image stretch saturating some of the ground features. To see more detail in these saturated ares, let's set the `Number of standard deviations` to `4`. Stretch standard deviation 4 Now you can see detail in some of the fields that were previously saturated. From the symbology pane, different stretch types and band combinations can be used to view Wyvern's hyperspectral data in different ways! --- [Source: https://knowledge.wyvern.space/docs/tutorials/envi/dragonette_data_in_envi] # Using Dragonette Imagery Data In ENVI Software ) ## Introduction The purpose of this tutorial is to provide a quick high-level introduction to some of the different ways in which Wyvern's satellite imagery data products acquired by their Dragonette satellites can be utilized in the ENVI® desktop software application while highlighting a handful of real-world analysis use cases. Rather than an exhaustive walkthrough, this tutorial covers some of the most popular basic functionality available in the ENVI software that is compatible with Wyvern's imagery products, including data ingest, preprocessing, interactive analysis, and automated algorithms for hyperspectral data exploitation. Finally, Wyvern does not explicitly endorse any specific commercial software tools and cannot provide technical assistance on the utilization of the ENVI software, if you have any specific questions on the ENVI software please contact your NV5 technical support representative. ### ENVI® Software Overview The ENVI® geospatial analysis software is developed, distributed, maintained, and supported by [NV5 Geospatial](https://www.nv5geospatialsoftware.com/). It is an industry standard image processing and geospatial data analysis application used by image analysts and data scientists around the world to extract timely, reliable, and accurate information from remote sensing imagery. The focus of this tutorial is on the modern ENVI v5+ desktop software interface and this tutorial does not cover the legacy ENVI Classic application. Figure 1 – Wyvern Dragonette-1 imagery data product loaded into the ENVI software (CIR) ENVI is a leader in spectral image processing and analysis, with top tools for multi and hyperspectral data including spectral target detection and identification. These tools are based on established, scientific methods for spectral analysis, using pixel responses at different wavelengths to obtain information about the materials within each pixel. They are used to detect objects, map features, measure vegetation health, identify minerals, detect marine debris, measure pollution, analyze wildlife habitats, map oil slicks, evaluate water quality, mitigate wildfires, detect methane leaks, perform supervised land use land cover (LULC) classification, and provide novel intelligence insights for numerous peace & security use case applications. ## Data Ingest & Display The first step when working with Dragonette imagery products is to ingest the raster dataset into the ENVI desktop software. The ENVI software has built-in smart data ingest readers which have intelligent support for a wide variety of [Raster Formats](https://www.nv5geospatialsoftware.com/docs/SupportedFormats.html#Raster) including Wyvern's Dragonette imagery data products. The optimal method for opening Wyvern imagery data products will depend on the version of ENVI desktop software that is being utilized. ### ENVI v6.0 (or older) In order to open the Wyvern imagery data product into the ENVI desktop software application use the `File` -> `Open…` menu option then select the imagery data product raster TIFF file or drag-n-drop this `*.tiff` file onto the ENVI display window. While the TIFF file format has band center wavelengths metadata embedded into the `*.tiff` file which are ingested by ENVI there is not currently a dedicated TIFF tag for storing wavelength units which are important to assign to the raster dataset in order to unlock certain downstream ENVI software functionality. Consequently, the recommended process of opening a Wyvern imagery dataset into ENVI v6.0 (or older) software is the following multi-step process: 1. `File` -> `Open…` menu option then select the imagery data product raster TIFF file (or drag-n-drop this `*.tiff` file onto the ENVI display window). 2. In the `Layer Manager` panel on left-hand side of ENVI application either double-click on the `wyvern_dragonette-*` image layer or right-click on the layer and in the resulting context menu select `View Metadata`. 3. At the bottom-right corner of the resulting `View Metadata` pop-up window press the `Edit Metadata` button. 4. At the upper-left corner of the resulting `Edit ENVI Header` pop-up window press the `+ Add…` button. 5. In the resulting `Add Metadata Items` pop-up window scroll down in the list and select `Wavelength Units` then press the `OK` button to dismiss the dialog. 6. Back in the `Edit ENVI Header` dialog navigate to the `Spectral` tab and make sure the `Wavelength Units` field drop-down list is set to `Nanometers` then press `OK` button with the `Display Result` box checked (on) which will reload the image into the display window. ### ENVI v6.1 (or newer) In order to open the Wyvern imagery data product into the ENVI desktop software application use the `File` -> `Open As` -> `Optical Sensors` -> `Wyvern` menu option then select the STAC metadata JSON file (in the same subfolder next to the TIFF image file) or drag-n-drop this `*.json` metadata file onto the ENVI display window. While using the menu option `File` -> `Open…` or drag-n-drop opening of the TIFF image file will still work in the ENVI v6.1 (or newer) application the dedicated Wyvern raster data ingest method using the JSON metadata file provides enhanced downstream ENVI software functionality. Consequently, the recommended process of opening a Wyvern imagery dataset into ENVI v6.1 (or newer) software is a much simpler single-step process: 1. `File` -> `Open As` -> `Optical Sensors` -> `Wyvern` menu option then select the STAC metadata JSON file (or drag-n-drop this `*.json` file onto the ENVI display window). 2. In the `Layer Manager` panel on left-hand side of ENVI application either double-click on the `wyvern_dragonette-*` image layer or right-click on the layer and in the resulting context menu select `View Metadata`. Figure 2 – Raster metadata information for a Wyvern Dragonette Level-1B imagery dataset ### Color Band Combinations By default when a Wyvern Dragonette imagery data product is opened into the ENVI software it will automatically display a True Color (RGB) band combination image rendering of the raster dataset in the main display window at 100% (1:1) resolution positioned at the center of the image scene. A detailed overview of how to use the ENVI software user interface and the display window navigation controls (Pan, Zoom, Rotate, Go To, etc.) is beyond the scope of this tutorial so for more information on the ENVI user interface basics please refer to NV5's online documentation for the ENVI software [`Display Tools`](https://www.nv5geospatialsoftware.com/docs/DisplayTools.html). In order to generate the RGB true color image display the ENVI software will automatically select the spectral bands with center wavelengths closest to the Red, Green, and Blue portions of the visible electromagnetic spectrum. Furthermore, the ENVI software has other [`Dynamic Band Selection`](https://www.nv5geospatialsoftware.com/docs/RasterLayers.html#BandSelect) controls that are easily customized via the panel in the lower-right hand corner of the main ENVI user interface (which is available at the bottom of the Layer tab whenever a raster is the selected layer in the Layer Manager). Since Dragonette imagery data has spectral bands that cover the visible and near-infrared (VNIR) wavelengths there are two predefined color band combinations that work well with Wyvern datasets: `True Color` (RGB) and `Color Infrared` (CIR). Furthermore, users can customize the image display band combination via click-and-drag of the individual blue & green & red slider bars in order to interactively select the bands for each color channel. 1. Click on the `wyvern_dragonette-*` image layer in the `Layer Manager` panel to make sure it is the selected active layer. 2. In the lower-right hand corner of the ENVI application window a `Layer` tab will be displayed that contains useful properties of the Dragonette imagery raster dataset. 3. At the bottom of this `Layer` tab panel click the dropdown arrow button to change the image display between `True Color` or `Color Infrared` predefined color band combinations in the resulting pop-up menu and observe how the image rendering is updated in the main display window. Experiment with other band combinations by moving the Blue & Green & Red tick marks on the color slider: Figure 3 – Comparison of True Color (RGB) and Color Infrared (CIR) band combinations ### Stretch Types & Stretch Extents By default, the ENVI software will automatically display imagery datasets with a default [Stretch Type](https://www.nv5geospatialsoftware.com/docs/DisplayTools.html#StretchTypes), and for Wyvern imagery with a 32-bit floating-point (float32) pixel data type the `Linear 2%` quick stretch type will be applied. While the Linear 2% stretch should provide a reasonable initial image visualization it is best suited to the mid-range pixels, since it maximizes contrast by saturating the lowest 2% of pixels to all black and the highest 2% to all white while stretching the mid-range pixel values in between. Consequently, in order to effectively visualize features in the darkest and brightest range of pixels while still maintaining visual separation of nuanced pixel radiance intensity differences end users may benefit from selection of alternative stretch types and/or the stretch on extent options available in the ENVI software. 1. Open a Wyvern Dragonette imagery data product and display a color band combination image into the view window. By default, the ENVI software will display the image with a stretch type of `Linear 2%` which will be the selection in the stretch type drop-down list: 2. On the main toolbar at the top of the ENVI software window click on the stretch type drop-down list and experiment by selecting other stretch types. In particular, the `Optimized Linear` stretch type will avoid the saturation of low and high pixel intensities at the tail ends of the image histogram thereby retaining visual fidelity in the darker and brighter regions of the imagery: 3. To the left of the stretch type drop-down list are three buttons that control the imagery area extent that is used to calculate the statistics for the current image stretch applied in the view display window. By default, the button will be selected for the `Stretch on Full Extent` mode which means the image stretch is based on statistics calculated on the data histogram of pixel values from the entire full extent of the imagery dataset. 4. In order to optimize the display stretch for the portion of the image currently displayed in the view window experiment with the other two stretch extent options: 1. Press the button to select `Stretch on View Extent` mode which updates the image stretch for the current view extent a single time 2. Press the button to select `Stretch on View Extent with Auto Update` mode which continuously updates the image stretch based on the current view extent while Pan & Zoom & Rotate movement occurs Figure 4 – Image with Optimized Linear stretch type (left) and stretch on view extent (right) ### Brightness & Contrast & Sharpen Once the desired image display stretch type and stretch extent has been selected the ENVI software also includes additional visual [Enhancement Tools](https://www.nv5geospatialsoftware.com/docs/DisplayTools.html#EnhancementTools) which can be applied by leveraging the brightness, contrast and sharpen sliders on the main toolbar at the top of the ENVI software window. All three of these image display enhancements can be used to improve the visual appearance of the imagery rendering within the active view window. 1. Start by selecting the `Linear` option from the stretch type drop-down list. 2. Locate the `Brightness` slider and drag the tick marker to a brightness value of `65`. This will slightly increase the overall brightness of the current image display. 3. Locate the `Contrast` slider and drag the tick marker to a contrast value of `0`. This will eliminate any extra contrast enhancement in order to display the original imagery. 4. Locate the `Sharpen` slider and drag the tick marker to a sharpen value of `100`. This will apply the maximum sharpening which will enhance the edges of smaller features within the imagery. Figure 5 – Image display with no Contrast (left) and increased Brightness & Sharpen (right) ## Statistical Analysis The ENVI software provides a wide variety of statistical analysis techniques for imagery raster datasets that can quickly & easily identify anomalous image pixels or segregate spectrally distinct regions. ### RX Anomaly Detection The ENVI software includes a simple [RX Anomaly Detection](https://www.nv5geospatialsoftware.com/docs/RXAnomalyDetection.html) tool that runs the Reed-Xiaoli Detector (RXD) algorithm which identifies pixels that are spectrally distinct from the overall imagery dataset background and thus represent anomalies within the image scene. In order for the RXD algorithm to be effective the anomalous targets must be relatively small in comparison to the overall spatial extent of the image scene. Results from the RXD analysis tool are unambiguous and have proven highly effective in detecting subtle spectral features within any given hyperspectral imagery dataset. The output raster generated by the RXD algorithm can also highlight anomalous scan lines but these do not affect the detection of other valid spectrally anomalous pixels within the image scene. 1. Within the Toolbox panel on the right-hand side of the ENVI software window type `rx` string in the blank search text box, press the `Enter` key on the keyboard, and double-click on the resulting `RX Anomaly Detection` tool that is highlighted. 2. Make sure to select the Wyvern Dragonette imagery dataset as the `Input Raster`, set the `Suppress Vegetation` parameter to the appropriate setting depending on whether vegetation should be considered anomalous by the RXD algorithm, then specify an output raster filename. 3. Press the `OK` button to run the RXD algorithm and the resulting output anomaly raster dataset will be generated then displayed as a single-band grayscale image with higher pixel values representing more spectrally anomalous pixels. 4. Within the Toolbox panel reset the search and enter `band roi` string in the blank search text box, press the `Enter` key, and double-click on the resulting `Band Threshold to ROI` tool that is highlighted. 5. In the resulting `Choose Threshold` dialog click-n-drag the maximum line on the anomaly raster histogram plot then move the maximum line all the way to the right-hand side of the plot window to set the `Max Value` for the region-of-interest (ROI) as the maximum RXD raster dataset pixel value. 6. Set the `Min Value` for the ROI by typing a minimum anomalous pixel intensity value of `20` in the text box then press the `Enter` key on the keyboard. 7. Once the band threshold ROI is generated in the `Data Manager` pop-up window right-click on the `New ROIs` item in the list and select `Load`. Figure 6 – Band threshold to ROI dialog with selection of high anomaly detection pixels The resulting ROI that was generated by thresholding the highest anomalous pixel values from the RXD algorithm output raster will be displayed as bright red color on top of the current active view image display. Figure 7 – RGB color image (left) and red anomaly pixels identified by RXD algorithm (right) ### 2D Scatter Plot The ENVI software includes a convenient [2D Scatter Plot](https://www.nv5geospatialsoftware.com/docs/using_envi_ScatterPlots.html) visualization tool that plots the pixel intensity values from one spectral band versus another spectral band which helps illustrate the degree of correlation between the two bands. It also provides a variety of interactive data exploration capabilities. For example, the scatter plot tool provides the ability to interactively view the distribution of pixel values within a patch moved under the cursor by clicking on the image display (called 'dancing pixels') along with region-of-interest drawing capabilities within the scatter plot in order to perform a simple supervised classification of the hyperspectral dataset. 1. With an imagery dataset displayed and selected in the `Layer Manager` panel select `Display` -> `2D Scatter Plot` from the main menu or press the button on the toolbar at the top of the main ENVI software user interface window. 2. Experiment with the interactive band selection sliders available next to the X (horizontal) and Y (vertical) axes in order to select different spectral bands for the scatter plot visualization. Note that if the same band is selected for both X and Y the resulting scatter plot will be a simple linear line with a uniform slope since the same band is perfectly correlated with itself. 3. Setup the scatter plot tool visualization with a selection of the green wavelength `Band 0.535 µm` for the X axis band and infrared wavelength `Band 0.799 µm` for the Y axis band. 4. From the `Scatter Plot Tool` dialog menu select `Options` -> `Patch Size` -> `25` in order to increase the pixel patch size for the interactive dancing pixels visualization. 5. Back on the main ENVI software window click-n-drag within the active view display window to move the pixel patch box around in the image then see the corresponding **dancing pixels** highlighted in bright red color within the 2D scatter plot window. 6. Within the `Scatter Plot Tool` dialog window click-and-drag to interactively draw a polygon around an isolated section of the scatter plot pixels in order to define region-of-interest (ROI) class. In order to classify unique distinct features in the imagery dataset draw ROI class polygons around the separate isolated clusters of pixels in the scatter plot. 7. Repeat the process of drawing separate distinct ROI polygons for multiple separate classes by pressing the button to create a new class then draw a corresponding polygon within the scatter plot and note that any given pixel in the scatter plot can only belong to one individual ROI class. Figure 8 – 2D scatter plots with dancing pixels (left) and interactive class definition (right) By creating multiple separate region-of-interest (ROI) classes in the 2D scatter plot dialog a simple supervised classification is effectively generated with each class color displayed on top of the source imagery dataset. These ROIs can also be used in subsequent supervised classification algorithms and processing workflows such as the `Classification Workflow` available in the toolbox panel on the right-hand side of the ENVI software window. Figure 9 – Multiple region-of-interest (ROI) classes defined using the 2D scatter plot tool ### ICA Transformation The ENVI software includes the [Independent Component Analysis (ICA)](https://www.nv5geospatialsoftware.com/docs/IndependentComponentsAnalysis.html) transformation algorithm which highlights the spectrally distinct features in the hyperspectral imagery data by amplifying the signal and suppressing the noise. The independent component analysis (ICA) algorithm decomposes data into independent signals and is a powerful processing technique for hyperspectral data in order to unmix the statistical independent endmembers. The ICA transformation can distinguish features of interest even when they occupy only a small portion of the pixels in the image. 1. Within the Toolbox panel on the right-hand side of the ENVI software window type `forward ica new` string in the blank search text box, press the `Enter` key on the keyboard, and double-click on the `Forward ICA Rotation New Statistics and Rotate` tool. 2. Select a Wyvern Dragonette hyperspectral imagery dataset in the resulting `Independent Components Input File` dialog then press the `OK` button. 3. Leave all ICA algorithm options configured with their default settings and specify an output ICA transformation raster dataset filename by pressing the `Choose` button to the right of the `Enter Output Filename` label. 4. Press the `OK` button to run the ICA algorithm and the resulting output ICA transformation raster dataset will be generated. 5. Display the ICA output raster dataset as RGB color using the lowest three bands which represent the most spectrally distinct and statistically independent data components of the input hyperspectral imagery dataset. Figure 10 – The first three bands of ICA transformation displayed as RGB color image (right) ## Atmospheric Correction Wyvern's Dragonette imagery data products delivered with Level-1B (L1B) processing level have been radiometrically corrected to pixel units of at-sensor top-of-atmosphere (TOA) radiance based on spacecraft location plus pointing along with solar conditions at time of data acquisition. Since the L1B imagery data product is delivered as a raster dataset with 32-bit floating-point (float32) data type where the pixels represent top-of-atmosphere (TOA) radiance the imagery data is already in units of `W / (m² * sr * µm)` with no need to apply any scaling factor. Consequently, the L1B imagery data product can be immediately visualized and directly analyzed using the statistical analysis techniques described in the previous sections of this tutorial. However, since the L1B imagery dataset's pixel values represent the solar illumination intensity as recorded by the hyperspectral imaging sensor on the Dragonette satellite in space the image data includes the impact of gases and aerosols in the atmosphere. In other words, the sunlight radiation has passed through the atmosphere twice before the final solar illumination intensity is recorded by the hyperspectral imaging sensor. Consequently, if the goal is analysis of the materials on the surface of the Earth a common pre-processing workflow is to first apply an atmospheric compensation correction to the L1B imagery product in order to remove the impact of atmospheric gases & aerosols thereby converting the image raster dataset pixels from top-of-atmosphere (TOA) radiance units to bottom-of-atmosphere (BOA) surface reflectance units. Fortunately, the ENVI software has a dedicated [Atmospheric Correction Module (ACM)](https://www.nv5geospatialsoftware.com/docs/AboutAtmosphericCorrectionModule.html) with several atmospheric compensation modeling tools for this conversion. The atmospheric correction algorithms available in the ENVI software model the amount of gases, water vapor, and aerosols in the atmosphere based on the input source hyperspectral imagery dataset and other known metadata parameters such as the geographic location, acquisition time, solar elevation angle, etc. The computed atmospheric properties are then used to constrain highly accurate models of radiation transfer in order to remove the effects of the atmosphere and generate a new output raster imagery dataset where the pixel values represent an accurate estimate of the true surface reflectance at the bottom-of-atmosphere (BOA). ### QUick Atmospheric Correction (QUAC) One of the most powerful atmospheric compensation modeling algorithms available in the ENVI software for pre-processing of Wyvern's Level-1B (L1B) imagery data product is the [QUick Atmospheric Correction (QUAC)](https://www.nv5geospatialsoftware.com/docs/quac.html) tool. The ENVI software QUAC tool provides atmospheric correction of Wyvern's hyperspectral imagery datasets with spectral coverage across the visible and near-infrared (VNIR) wavelengths. QUAC determines the atmospheric correction parameters directly from the observed pixel spectra in an image scene without the need for the end-user to manually input additional ancillary information. QUAC is based on the empirical finding that the average reflectance of diverse material spectra is not dependent on each scene so processing is much faster compared to first-principles atmospheric modeling methods. QUAC allows for any imaging off-nadir angle or solar elevation angle and QUAC works best when there are multiple diverse materials in an image scene with sufficiently dark pixels to allow for a good estimation of the baseline spectrum. 1. Within the Toolbox panel on the right-hand side of the ENVI software window type `quac` string in the blank search text box, press the `Enter` key on the keyboard, and double-click on the `QUAC - Quick Atmospheric Correction` tool. 2. Select a Wyvern Dragonette imagery dataset and press the `OK` button. 3. In the resulting `QUAC - Quick Atmospheric Correction` dialog window make sure the `Sensor Type` drop-down list is set to `Generic / Unknown Sensor` which is the appropriate selection for Wyvern Dragonette imagery datasets. 4. Specify an `Output Raster` dataset filename and press the `OK` button to run QUAC. Figure 11 – Configuration of the QUAC - Quick Atmospheric Correction processing tool ## Spectral Analysis The ENVI software provides a wide variety of simple spectral analysis techniques for exploitation of hyperspectral imagery datasets. For the most accurate analysis of the materials on the surface of the Earth, the Wyvern L1B imagery data product should first be atmospherically corrected as described in the previous Atmospheric Correction section. ### Spectral Profiles The ENVI software provides a wide variety of interactive data plotting visualization tools where a common way to analyze hyperspectral imagery datasets is to extract spectral profiles (also known as Z profiles since they plot data values in the Z spectral dimension of a raster hypercube dataset). The [Spectral Profile (Z Profile)](https://www.nv5geospatialsoftware.com/docs/Profiles.html#SpecProfile) tool in the ENVI software plots the pixel intensity spectrum for all spectral bands for any given pixel selected by the end-user clicking on selected image layers within the current active view window. ENVI uses the imagery product metadata and any custom header information to automatically scale and label the spectral profile plot. If multiple layers are displayed in the active view ENVI plots a Spectral Profile from each layer, to turn off the plot for an individual layer simply disable the check box for that layer in the Layer Manager. 1. With an imagery dataset displayed and selected in the `Layer Manager` panel select `Display` -> `Profiles` -> `Spectral` from the main menu or press the button on the toolbar at the top of the main ENVI software user interface window. 2. Make sure the ENVI software navigation controls are in active `Select` mode by selecting the button on the top toolbar or press the `F5` key on the keyboard. 3. Using the mouse **single-click** on different pixels within the image display view window and the corresponding spectrum will be plotted as a line within the Spectral Profile dialog. The resulting spectral profile line plot represents the intensity of raster data values for the selected pixel throughout all of the spectral bands in the hyperspectral imagery data cube. 4. In order to analyze how the spectral signatures vary within the image scene use the mouse to **click-n-drag** across the image display view window then watch how the spectral profile line plot dynamically updates. 5. In order to collect multiple spectral profiles for different pixels within the image hold down the **Shift key** on the keyboard then use the mouse to **left-click repeatedly** within the image display view window on different areas of interest within the image scene and the corresponding spectral profile will be plotted for each pixel clicked. 6. On the right-hand side of the `Spectral Profile` dialog press the `>` bar button to expand the plot properties panel in order to see the list of collected spectra for the multiple pixels and select the `Plot Stats` tab to see useful data statistics for each pixel spectrum. Figure 12 – Multiple spectral profiles generated with useful statistics for each pixel spectrum The spectral profile plotting tool also provides a convenient way to visualize and inspect raster data pixel values from multiple imagery datasets in order to compare their respective spectral signatures. For example, an original Wyvern Level-1B (L1B) imagery dataset with top-of-atmosphere (TOA) radiance pixel units can be compared with the atmospherically compensated output result from the QUick Atmospheric Correction (QUAC) tool described in the previous section which generates a new imagery dataset with bottom-of-atmosphere (BOA) surface reflectance pixel units. 1. From the main ENVI software menu select `Views` -> `Two Vertical Views` or press `Ctrl+2` on the keyboard. 2. Using the instructions in the ENVI software documentation for [`Multiple Views`](https://www.nv5geospatialsoftware.com/docs/views.html#multiview) load the original Wyvern L1B imagery dataset into the first (left) view window then load the QUAC output imagery dataset into the second (right) view window. 3. From the main ENVI software menu select `Views` -> `Link Views` then in the resulting `Link Views` dialog select the `Geo Link` radio button option, press the `Link All` button, then press `OK` button to link the views and dismiss the dialog. This will ensure the two views are geographically linked so that each view displays the same location and extent of imagery. 4. Within the `Layer Manager` panel on left-hand side of ENVI software window select the original Wyvern L1B imagery dataset raster layer then press the button on the main toolbar in order to launch a new `Spectral Profile` plot for this dataset. 5. Within the `Layer Manager` panel on left-hand side of ENVI software window select the QUAC output result imagery dataset raster layer then press the button on the main toolbar in order to launch a new `Spectral Profile` plot for this dataset. 6. Make sure the ENVI software navigation controls are in active `Select` mode by selecting the button on the top toolbar or press the `F5` key on the keyboard. 7. To compare spectral signatures from both imagery datasets use the mouse then **click-n-drag** across one of the image display view windows then watch how the spectral profile line plots dynamically update. This spectral profile comparison also illustrates how the Level-1B imagery dataset's slightly lower pixel intensity values in the 764 nm center wavelength band due to sunlight transmission impedance caused by oxygen absorption in the atmosphere has been corrected in the QUAC tool processing result output raster imagery dataset with surface reflectance units. Figure 13 – Spectral profile comparison of original Wyvern Level-1B TOA radiance data (left) and the QUAC atmospheric correction output with BOA surface reflectance data (right) ### Material Identification The ENVI software includes a wide variety of advanced spectral analysis tools that can be used to extract useful insights from Wyvern's Dragonette hyperspectral imagery data products. One such capability is the [Material Identification Tool](https://www.nv5geospatialsoftware.com/docs/materialidentificationtool.html), which compares any given imagery dataset pixel's unknown spectral profile against the spectra of known materials in ENVI's rich collection of reference [Spectral Libraries](https://www.nv5geospatialsoftware.com/docs/spectraltoolsreferences.html#Spectral), then ranks the similarity between the two using common spectral similarity algorithms. The tool can be found in the expanded plot properties panel on the right-hand side of any given Spectral Profile plot dialog window. 1. Open and display the **QUAC atmospheric correction** output processing result raster dataset so the imagery data being analyzed has bottom-of-atmosphere (BOA) surface reflectance pixel units. 2. Launch a new Spectral Profile for the QUAC-processed imagery dataset by selecting `Display` -> `Profiles` -> `Spectral` from the main menu or press the button on the toolbar at the top of the main ENVI software user interface window. 3. On the right-hand side of the `Spectral Profile` dialog press the `>` right arrow bar to expand the plot properties panel then press the button labeled `Identify` in order to launch the `Material Identification` tool. 4. Press the `Select Library` drop-down list button and select an appropriate spectral library `*.sli` file based on the potential materials contained in the imagery dataset. For example, the `usgs_v7` -> `vegetation_aviris.sli` spectral library contains signatures for common vegetation materials derived from AVIRIS hyperspectral data which covers the same visible and near-infrared (VNIR) wavelengths as Wyvern Dragonette imagery datasets so it makes an excellent reference library for vegetative material identification. 5. Using the mouse **single-click** on the image display window to select a specific target pixel of interest then note how the `Material Identification` tool displays a top-down sorted table of potential spectral matches from the selected reference library. Figure 14 – Material identification based on USGS v7 spectral library showing high likelihood that selected imagery dataset pixel matches spectral signatures for Lodgepole Pine tree. ### Spectral Indices Since Wyvern's hyperspectral imagery datasets contain dozens of contiguous spectral bands with narrow bandwidths across visible and near-infrared (VNIR) wavelengths a very common and extremely powerful analysis technique is the extraction of [Spectral Indices](https://www.nv5geospatialsoftware.com/docs/spectralindices.html) by applying specific predefined band algebra equations. The remote sensing science community has defined numerous spectral indices which leverage specific algebraic equations for band math arithmetic in order to highlight different features and properties within an Earth observation imagery dataset. Spectral indices are combinations of spectral reflectance from two or more wavelengths that indicate the relative abundance of features of interest such as fire burn severity, geologic minerals, man-made objects, built-up features, water & snow, soil & mud, and a wide variety of advanced vegetation properties. The ENVI software includes a convenient analytics tool for generating almost a hundred different spectral indices as delineated in the [Alphabetical List of Spectral Indices](https://www.nv5geospatialsoftware.com/docs/alphabeticallistspectralindices.html) including several specialized [Narrowband Greenness](https://www.nv5geospatialsoftware.com/docs/narrowbandgreenness.html) vegetation indices which are specifically designed for advanced vegetative analysis of hyperspectral imagery datasets. 1. Open and display the **QUAC atmospheric correction** output processing result raster dataset so the imagery data being analyzed has bottom-of-atmosphere (BOA) surface reflectance pixel units. 2. Within the Toolbox panel on the right-hand side of the ENVI software window type `indices` string in the blank search text box, press the `Enter` key on the keyboard, and double-click on the `Spectral Indices` tool. 3. Select the QUAC output raster dataset and press `OK` button. 4. In the resulting `Spectral Indices` dialog window select the desired indices in the `Index` list in order to generate the new index raster(s) from the atmospherically corrected imagery dataset: 5. Once the spectral indices processing task is finished running open the `Data Manager` dialog, right-click on a spectral index band in the resulting output raster dataset, and select `Load Grayscale`. 6. Within the `Layer Manager` panel on left-hand side of ENVI software right-click on the spectral index raster layer, select `Change Color Table` -> `More…`, then within the `Change Color Table` dialog press the button to select a different color table. Figure 15 – Display of several spectral indices with different color tables applied ## Advanced Analysis ### Target Detection Since the ENVI software includes a rich collection of several spectral libraries that contain the spectral signatures of a wide variety of natural and man-made materials these spectral libraries can be used as a baseline reference to perform spectral target detection analysis. The ENVI software includes a convenient guided multi-step [Target Detection Workflow](https://www.nv5geospatialsoftware.com/docs/TargetDetectionWorkflow.html) that locates regions and objects within a hyperspectral imagery dataset whose pixels are spectrally similar to target spectral signatures from a reference spectral library. In this tutorial a spectral signature for light gray concrete road material will be used from a reference spectral library included with the ENVI software. 1. From the main ENVI software menu select `Display` -> `Spectral Library Viewer`. 2. In the blank `Search` text box enter the string `concrete` and press the `Enter` key on the keyboard. 3. Navigate into the `usgs_v7` subfolder, select the `artificial_asdfr.sli` spectral library file for artificial materials. 4. Click on the `Concrete_GDS375_Lt_Gry_Road_ASDFRa` item in the list and the corresponding spectral signature for this concrete material will be plotted. Figure 16 – Spectral signature for light gray road concrete material from USGS v7 library Once an appropriate reference spectral signature has been identified for the material of interest the workflow can be executed against an atmospherically corrected hyperspectral imagery dataset. 1. Open and display the **QUAC atmospheric correction** output processing result raster dataset so the imagery data being analyzed has bottom-of-atmosphere (BOA) surface reflectance pixel units. 2. Within the Toolbox panel on the right-hand side of the ENVI software window type `target detect` string in the blank search text box, press the `Enter` key on the keyboard, and double-click on the `Target Detection Workflow` tool. 3. In the workflow `Select Data` step to the right of the blank `Input Raster` text box press the `Browse…` button, select the QUAC output raster dataset and press `OK` button, then press the `Next >>` button to proceed to the next step in the workflow. 4. In the workflow `Select Signatures` step make sure the `Target` subfolder list item is selected then press the `Import` button and from the drop-down list menu select `From Spectral Library…`. 5. In the resulting `Import from Spectral Library` dialog within the top `Library` section list scroll down and select the `usgs_v7\artificial_asdfr.sli` spectral library file. Next, in the blank `Search…` text box to the right of `Spectra` label enter the text string `concrete` then select the `Concrete_GDS375_Lt_Gry_Road_ASDFRa` item in the list, press `OK` button to dismiss the pop-up dialog, then press the `Next >>` button to proceed to the next step in the workflow. 6. In the workflow `Image Transform for Dimensionality Reduction` step select `Skip this step` radio button then press the `Next >>` button to proceed to the next step in the workflow. 7. In the workflow `Target Detection` step click the `Method` drop-down list, select `Spectral Angle Mapper Classification` algorithm from the menu, then press the `Next >>` button to proceed to the next step in the workflow. 8. In the workflow `Threshold` step interactively drag the red threshold bar on the SAM classification histogram to an appropriate level aligned with peak spectral target detections then press the `Next >>` button to proceed to the next step in the workflow. 9. In the workflow `Smooth` step leave all parameters set to their default settings then press the `Next >>` button to proceed to the next step in the workflow. 10. In the workflow `Export Final Result` step specify output filenames for both the `Export Classification Raster` and `Export Shapefile` then press the `Finish` button to complete execution of the target detection workflow. Figure 17 – CIR image (left) and spectral target detection of light gray road concrete (right) Once the spectral target detection workflow is finished running the output classification raster and vector shapefile datasets can be displayed to illustrate the regions in the imagery where there is a high likelihood match to the spectral signature for light gray road concrete material from the reference spectral library. ### Supervised Classification The ENVI software includes a wide variety of advanced [Classification Tools](https://www.nv5geospatialsoftware.com/docs/ClassificationTools.html) that can be used to perform a supervised classification of hyperspectral imagery datasets in order to generate a map of different land use land cover (LULC) types. The supervised classification process in the ENVI software is a multi-step workflow which involves the manual creation of supervised training data for the different class types, execution of a supervised classification algorithm, post-classification cleanup, and raster-to-vector conversion for processing result output to common GIS mapping formats such as Shapefile. The ENVI software provides a powerful [Region of Interest (ROI) Tool](https://www.nv5geospatialsoftware.com/docs/RegionOfInterestTool.html) for definition of supervised training data along with a convenient guided [Classification Workflow](https://www.nv5geospatialsoftware.com/docs/classification.html) that walks users through these multiple steps in order to generate GIS mapping layers. 1. Open and display the **QUAC atmospheric correction** output processing result raster dataset so the imagery data being analyzed has bottom-of-atmosphere (BOA) surface reflectance pixel units. 2. While the QUAC-processed imagery dataset is selected in the `Layer Manager` panel on the left-hand side of the ENVI software window launch the `ROI Tool` by selecting `File` -> `New` -> `Region of Interest` from the main menu or press the button on the toolbar at the top of the main ENVI software user interface window. 3. Within the `Region of Interest (ROI) Tool` pop-up dialog window press the button to initialize a new ROI and in the text box to the right of the `ROI Name` label change the name of the *ROI #1* to a more suitable name for the training data class type being defined (e.g., *Open Water*). 4. Follow the instructions in the ENVI software documentation on how to utilize the [`Region of Interest (ROI) Tool`](https://www.nv5geospatialsoftware.com/docs/RegionOfInterestTool.html) in order to define separate ROI classes for each unique land use land cover (LULC) type in the geographic region where the atmospherically corrected hyperspectral imagery dataset is located. 5. One of the simplest ways to quickly and effectively create supervised training data ROIs for remote sensing imagery is to leverage the [Create ROIs from Pixels](https://www.nv5geospatialsoftware.com/docs/RegionOfInterestTool.html#ROIPixel) functionality available when the `Region of Interest (ROI) Tool` dialog window is selected. 6. Repeat the process of creating a new ROI for every land use land cover (LULC) class type that should be classified using the supervised image classification workflow. 7. Once finished creating all of the desired ROIs save the training data to an XML format file on disk by selecting `File` -> `Save As…` from the `Region of Interest (ROI) Tool` dialog menu. Figure 18 – The ENVI software Region of Interest (ROI) Tool being used to manually create separate supervised training data ROIs for distinct land use land cover (LULC) class types Now that the [ROI Tool](https://www.nv5geospatialsoftware.com/docs/RegionOfInterestTool.html) has been used to generate supervised training data, the next step is to run the [Classification Workflow](https://www.nv5geospatialsoftware.com/docs/classification.html). 1. Within the Toolbox panel on the right-hand side of the ENVI software window type `classification workflow` string in the blank search text box, press the `Enter` key on the keyboard, and double-click on the `Classification Workflow` tool. 2. In the workflow `File Selection` step make sure the QUAC output raster dataset is selected for the `Input Raster File:` field then press the `Next >` button to proceed to the next step in the workflow. 3. In the workflow `Classification Type` step select the `Use Training Data` radio button then press the `Next >` button to proceed to the next step in the workflow. 4. In the workflow `Supervised Classification` step the previously defined ROIs should be automatically loaded into the `Training Data` list panel (if not then press the `Load Training Data Set` button and load the ROIs previously saved to the XML format file on disk). 5. Select the `Algorithm` tab, click the drop-down list at the top of the tab panel, then select `Spectral Angle Mapper` algorithm. 6. In the lower left-hand corner of the `Classification` workflow dialog check the `Preview` checkbox then use the main image display navigation controls to move the underlying image around to view a preview of the SAM image classification result within the preview inlay box: 7. If necessary make advanced adjustments to the `Maximum Spectral Angle` by either specifying a `Single Value` threshold or set `Multiple Values` individually for each ROI class. 8. Once satisfied with the SAM image classification configurations press the `Next >` button to proceed to the next step in the workflow. 9. In the workflow `Cleanup` step leave all parameters in their default settings and press the `Next >` button to proceed to the next step in the workflow. 10. In the workflow `Cleanup` step specify output filenames for both the `Export Classification Image` and `Export Classification Vectors` then optionally export the classification statistics to a text file on disk. 11. Once finished specifying the output filenames press the `Finish` button to complete execution of the classification workflow. Once the Classification Workflow processing steps have finished running the output classification raster and vectors will be automatically loaded into the active view display window where the land use land cover (LULC) classification results can be visualized. Figure 19 – CIR image (left) and supervised classification workflow processing results (right) --- [Source: https://knowledge.wyvern.space/docs/tutorials/python/loading_data] # Loading & Viewing Wyvern Data in Python Output color infrared image from the Python notebook Try it yourself Runs free in your browser with no install or setup, and pulls the imagery straight from our Open Data Program. Or view the notebook source on GitHub Open in Colab > **Note:** There are many ways to visualize raster data using Python. The combination of Rasterio + Matplotlib was chosen due to the simplicity & ability to build off of this code. ## Getting started - **Google Colab (easiest):** open the notebook using the "Open in Colab" badge above, then run the cells top to bottom. - **Local Jupyter:** open the notebook in JupyterLab or VS Code (Python 3.10+) and run it. Either way, the first cell installs the packages that are not already present: ```python %pip install -q numpy matplotlib requests rasterio ``` ## Image we're using Check out the image on the [Wyvern Open Data Catalog](https://opendata.wyvern.space/#/year/2024/wyvern_dragonette-001_20240930T070744_08fd7f5a_l2a/wyvern_dragonette-001_20240930T070744_08fd7f5a_l2a.json)! This is L2A surface reflectance data, meaning atmospheric effects have already been corrected for and pixel values represent an estimate of actual ground reflectance rather than raw at-sensor radiance. ## Band Combinations We'll be exploring the following band combinations! | Name | Red Channel | Green Channel | Blue Channel | | ---- | ----------- | ------------- | ------------ | | Natural Color (RGB) | 650nm | 550nm | 503nm | | Color Infrared (CIR) | 750nm | 650nm | 550nm | ## Code First things first, let's import all of our packages. ```python import json import numpy as np import rasterio import matplotlib.pyplot as plt ``` ## Download the Tutorial Files Let's download the GeoTIFF and its STAC JSON metadata! Run the cell below to download these files. Once the download is complete, they'll appear in the same folder as your notebook! ```python import requests from pathlib import Path IMAGE_URL = "https://wyvern-data.com/wyvern_dragonette-001_20240930T070744_08fd7f5a_l2a/wyvern_dragonette-001_20240930T070744_08fd7f5a_l2a.tiff" JSON_URL = "https://wyvern-odp.com/wyvern_dragonette-001_20240930T070744_08fd7f5a_l2a/wyvern_dragonette-001_20240930T070744_08fd7f5a_l2a.json" IMAGE_PATH = Path("wyvern_dragonette-001_20240930T070744_08fd7f5a_l2a.tiff") JSON_PATH = Path("wyvern_dragonette-001_20240930T070744_08fd7f5a_l2a.json") downloads = [ (IMAGE_URL, IMAGE_PATH), (JSON_URL, JSON_PATH), ] for url, path in downloads: if path.exists(): print(f"{path.name} File already here! Skipping download.") else: print(f"Downloading {path.name}...") response = requests.get(url, headers={"User-Agent": "Mozilla/5.0"}) response.raise_for_status() path.write_bytes(response.content) print("Downloads complete!") ``` ### Loading the Data We'll start by reading the scale factor from the metadata JSON file we downloaded earlier. Wyvern's L2A products represent surface reflectance, but the values are stored as scaled 16-bit integers to reduce file size. ```python # Open the JSON metadata file and extract the scale factor with open(JSON_PATH) as json_file: metadata = json.load(json_file) # The scale factor converts stored uint16 reflectance values into their floating-point representation scale_factor = metadata["assets"]["Cloud optimized GeoTiff"]["raster:bands"][0]["scale"] # if you don't have the JSON file, you can hardcode the scale factor: # scale_factor = 0.0001 print(f"Scale factor: {scale_factor}") ``` Now let's open the GeoTIFF using [rasterio](https://rasterio.readthedocs.io/en/stable/api/index.html). We'll replace any nodata pixels with `np.nan` so they are excluded from our calculations, then apply the scale factor to recover the floating-point reflectance values. ```python with rasterio.open(IMAGE_PATH) as image_file: # Convert to float so we can use np.nan and decimal values image_arr = image_file.read().astype(np.float32) nodata = image_file.nodata band_descriptions = image_file.descriptions # Replace nodata pixels with NaN so they are excluded from calculations if nodata is not None: image_arr[image_arr == nodata] = np.nan # Apply the scale factor to recover floating-point reflectance values image_arr *= scale_factor print(f"Image shape: {image_arr.shape}") print(f"Image bands: {', '.join(band_descriptions)}") ``` ## Building the Wavelength Map Unlike a standard RGB image, each band in a hyperspectral GeoTIFF corresponds to a specific wavelength. We can find this information in the `eo:bands` field under the GeoTIFF asset in the STAC JSON we downloaded earlier, which lists each band's center wavelength. We'll read these out and build a map from wavelength to band index so we can select bands by wavelength rather than hardcoding band numbers. ```python with open(JSON_PATH) as json_file: metadata = json.load(json_file) wavelength_map = {} for i, band in enumerate(metadata["assets"]["Cloud optimized GeoTiff"]["eo:bands"]): wavelength = round(band["center_wavelength"] * 1000) # convert from μm to nm wavelength_map[wavelength] = i print(wavelength_map) ``` ### Selecting Bands by Wavelength Now we'll write a small helper function that finds the nearest available band to any wavelength we ask for. This means if we request 550nm and the closest available band is 549nm, it'll find it automatically. ```python def get_band(wavelength_nm: int) -> np.ndarray: """ Returns the image band closest to the requested wavelength. Args: wavelength_nm: Target wavelength in nanometres. Returns: A 2D array containing reflectance values for the closest available wavelength band. """ available_wavelengths = list(wavelength_map.keys()) closest_wavelength = min(available_wavelengths, key=lambda w: abs(w - wavelength_nm)) print(f"Requested {wavelength_nm}nm -- using {closest_wavelength}nm") return image_arr[wavelength_map[closest_wavelength]] ``` Now that we have the image array loaded, we're good to plot it using [matplotlib](https://matplotlib.org/)! ```python # Select our R, G, and B bands by wavelength, then stack them into a [rows, columns, channels] # array which matplotlib expects rgb_arr = np.dstack([get_band(650), get_band(550), get_band(503)]).astype(np.float32) brightness_factor = 3 # Bring up the brightness a touch so we have a better visual rgb_arr_bright = rgb_arr * brightness_factor rgb_arr_bright = np.clip(rgb_arr_bright, 0, 1) # keep values within 0-1 # We specify a larger figsize to get a better resolution out of the image. You can # increase this further to further improve resolution! (at the expense of file size) fig, ax = plt.subplots(1, 1, figsize=(15, 10)) ax.imshow(rgb_arr_bright) plt.tight_layout() plt.show() ``` Output true-color image from the Python notebook We can also render the image in color infrared by selecting a different set of bands from the image. ```python # Now, let's do the same with color infrared! cir_arr = np.dstack([get_band(750), get_band(650), get_band(550)]).astype(np.float32) cir_arr_bright = cir_arr * brightness_factor cir_arr_bright = np.clip(cir_arr_bright, 0, 1) # keep values within 0-1 fig, ax = plt.subplots(1, 1, figsize=(15, 10)) ax.imshow(cir_arr_bright) plt.tight_layout() plt.show() ``` Output color infrared image from the Python notebook Well done! You now know how to download and view Wyvern images in Python! Jump to the next tutorial to make some super cool plots about hyperspectral pixels! --- [Source: https://knowledge.wyvern.space/docs/tutorials/qgis/loading_data] # Loading Data in QGIS & EnMAP Box ## Before you start... Before starting this tutorial, please ensure you have done the following: - QGIS has been downloaded installed on your machine. Visit [qgis.org](https://www.qgis.org/) to download QGIS. - You have [installed the EnMAP Box plugin into QGIS](https://enmap-box.readthedocs.io/en/latest/usr_section/usr_installation.html) - You have a Wyvern image you can load into QGIS. If you haven't downloaded one already, we'll be [using this image](https://opendata.wyvern.space/#/year/2024/wyvern_dragonette-001_20240930T070744_08fd7f5a/wyvern_dragonette-001_20240930T070744_08fd7f5a.json) for the next steps. :::info **Note:** Apple silicon devices may have trouble installing QGIS. For best results, use a Windows or Linux (x86) machine, or follow the instructions on the QGIS website to install the OSX version. ::: ## Loading into QGIS To load an image into QGIS, you have a number of different options: - Using the QGIS Browser - Using the QGIS Data Source Manager - Clicking the `Open Data Source Manager` button (`Ctrl + L` hotkey) - From the menu bar, accessing the `Add Raster Layer` option: `Layer` -> `Add Layer` -> `Add Raster Layer`, which brings up the Data Source Manager on the correct tab. (`CTRL + Shift + R` hotkey) - Dragging & dropping the image onto QGIS For this tutorial, we'll use the `Data Source Manager` option since it can be easily shown here. Go ahead and open up the Data Source Manager by clicking the `Open Data Source Manager` button below. Open Data Source Manager button Once you have opened up the Data Source Manager, it should look like the image below. Ensure you have clicked the `Raster` tab. Go ahead and click on the three dots (`...`) to bring up your file manager. QGIS Data Source Manager window Then, select a Wyvern geotiff to load into QGIS! Compatible images should show up when browsing folders. Selecting a Wyvern geotiff Some additional options have now appeared in the Data Source Manager. You can leave them all as-is, and click on the `Add` button to add the image to the current QGIS project. Selected image in Data Source Manager Congrats! You have successfully loaded an image into QGIS, you superstar, you! QGIS should look similar to the image below, which is rendering a Wyvern hyperspectral geotiff. :::info **NOTE:** QGIS does not know the correct red, green, and blue bands to display for our geotiff. This is why the image has washed-out colors. We'll fix that in the next section of the tutorial! ::: Wyvern geotiff rendered in QGIS ## Adjusting Band Rendering QGIS has no logic to select optimal bands to render an image. This means we'll have to do that manually! Luckily it is super easy to adjust the bands used to render the image. Below is a list of recommended bands for visualizations! *Note: generally, adjacent/close wavelengths to the below values will provide similar results.* | Name | Description | Red Band | Green Band | Blue Band | | ---- | ----------- | -------- | ---------- | --------- | | Natural Color (RGB) | Normal visualization | 650nm | 550nm | 503nm/490nm | | Color Infrared | Vegetation analysis | 750nm | 650nm | 550nm | ### Band Combinations To adjust our visualization settings, we have to open our `Layer Properties` window for our image/layer. To do so, right-click on our layer in the `Layers` panel, and select the `Properties` option, or double-click on the layer in the `Layers` panel. Selecting Layer Properties Now that we have the `Layer Properties` window open, ensure you have the `Symbology` tab selected. You should then be able to see dropdown boxes for red, green, and blue image channels. We can select the relevant band wavelengths from the table above, and select the closest band wavelengths in the image. **Note:** Each band is labeled `Band_XXX`, with `XXX` representing the centre wavelength of the band in nanometres. Symbology tab with band options In the image below, the relevant bands for a natural color (RGB) image have been selected. RGB bands selected for visualization Once you have selected your bands, click the `Apply` button to apply your changes and render your image! You should now see an RGB image displayed on the map! :::info At this point, the image may still look a little weird! This is okay, and due to the default contrast settings in QGIS! We'll adjust these values next! :::: RGB natural color rendering of image ### Adjusting Contrast Every image is slightly different in its lighting conditions and distribution of values within each band. QGIS has default contrast settings that work alright for most images, but definitely benefit from careful adjustment. With our image of irrigated fields in Saudi Arabia, we'll be turning our contrast adjustment to Min/Max to maximize image quality. To do this, first expand the `Min / Max Value Settings` option in the `Layer Properties -> Symbology` tab. You can then select the `Min / max` option for this image. ?> The `Mean +/- standard deviation` setting is a great option as well! Make sure to explore different options depending on the lighting conditions and ground cover of your image. Min/Max contrast selected Once you have selected the contrast option, hit the `Apply` button again to re-render the image. It should look fantastic! Wyvern hyperspectral image with correct contrast settings This image also looks great in color infrared! Wyvern hyperspectral image in color-infrared with correct contrast settings Well done! You now know how to load a Wyvern hyperspectral image into QGIS for analysis, and adjust rendering settings to render a natural color or color-infrared image! ## Loading into EnMAP Box Now that we have the image loaded into QGIS, let's do the same with EnMAP Box! ?> **Note:** QGIS and EnMAP-Box are separate applications. You do **not** need to load an image into QGIS before loading the same image into EnMAP Box. Before we can do that, we need to open EnMAP Box by selecting `Raster -> EnMAP-Box -> EnMAP-Box` from the menu bar ![EnMAP Box can be opened by selecting `Raster -> EnMAP-Box -> EnMAP-Box`](images/enmap_selection.jpg 'EnMAP-Box location in the menu bar :size=400') EnMAP-Box will then open as an additional window. At this point the window should be pretty empty. ![Inital EnMAP-Box window](images/enmap_empty.jpg 'Inital EnMAP-Box window :size=800') To load our image into EnMAP-Box, we can click the ![EnMAP Box Add Data Source](images/enmap_load_icon.jpg) EnMAP-Box load icon, which is located in the top-right of the window. This will open a file select dialog which will allow you to select our Wyvern image we downloaded in the previous steps. Our image should then appear in the `Data Sources` panel! We can right-click on the image and select `Open in new map -> True Color`, and EnMAP-Box will parse our image's metadata to select the most relevant band centres. How cool! ![Opening our raster in a new map](images/enmap_data_sources.jpg 'Opening our raster in True Color :size=800') Your image should look similar to the image below. Note the band centres are similar to the band centres we manually selected during QGIS import! ![Wyvern image loaded into EnMAP-Box](images/enmap_wyvern_image.jpg 'Wyvern image loaded into EnMAP-Box :size=1000') ## Adjusting Band Rendering If we wanted to change which bands we are viewing, we can **double-click** the image in the `Data Views` panel to open our normal QGIS `Symbology` window. EnMAP Box also provides an extremely useful `Color Space Explorer` tool which provides pre-defined band selections for specific analyses. To access the `Color Space Explorer` tool, navigate to the tool via the menu bar using `Tools -> Raster Visualizations -> Color Space Explorer`. ![Navigating to the Color Space Explorer](images/enmap_color_space_explorer.jpg 'Navigating to the Color Space Explorer :size=600') To use Color Space Explorer, first select our image from the top drop-down. Bands can then be selected in the `Band Rendering` tab, or you can just click the button you want from the bottom list of common band combinations! Wyvern's hyperspectral data is only Visible to Near-Infrared (VNIR), so only some of the combinations will work. These include: `Natural Color`, `False Color`, `Color Infrared`, and `Bathymetric`. ![Color Space Explorer window](images/enmap_color_space_explorer_window.jpg 'Color Space Explorer window :size=600') You're now able to load imagery into QGIS and EnMAP-Box! You rock! ![Color Space Explorer window](images/enmap_cir_render.jpg 'Wyvern image rendered in EnMAP-Box :size=1000') --- [Source: https://knowledge.wyvern.space/docs/tutorials/] # Tutorials and Code Wyvern is dedicated to making our imagery easy to use in your software of choice. This section of the Knowledge Centre is focused on easy-to-follow tutorials and examples. Each tutorial uses freely available imagery from Wyvern's [Open Data Program](https://opendata.wyvern.space/#/?.language=en). :::info Have questions, suggestions for new tutorials or GIS software, or just want to connect with us? Join us on our [Github Discussions](https://github.com/Nrevyw/wyvern-public-resources/discussions), or on [LinkedIn](https://www.linkedin.com/company/wyvern-space). We would love to chat more with you! ::: --- [Source: https://knowledge.wyvern.space/docs/documentation/getting_started/hyperspectral_basics] # Hyperspectral basics A hyperspectral image records light in many narrow, adjacent bands, so every pixel carries a spectrum rather than three colour values. The same technique is used in food sorting, pharmaceuticals, and astronomy. Here we look at it from orbit. This page starts with the light itself, then compares the three kinds of optical Earth observation imagery (panchromatic, multispectral, and hyperspectral) and what each can and can't tell you. ## Electromagnetic spectrum Light is electromagnetic radiation, and the electromagnetic spectrum is the full range of it, ordered by wavelength from gamma rays to radio waves. Optical Earth observation works in a narrow slice near the middle: visible light from roughly 400 to 700 nm, the near-infrared (NIR) from 700 to about 1,000 nm, and the shortwave infrared (SWIR) from about 1,000 to 2,500 nm. Visible and near-infrared together are called the VNIR. Sunlight reaches the surface, and each material absorbs some wavelengths and reflects others. Which ones depends on what the material is made of. A satellite measures the reflected light, and the pattern of what's missing is the information. The electromagnetic spectrum, with the visible range expanded. Optical satellites work in and just beyond the visible. `[4]` Wikipedia has a [more detailed explanation of the electromagnetic spectrum](https://en.wikipedia.org/wiki/Electromagnetic_spectrum) if you want the physics. ## Imagery types Optical sensors differ mainly in how they divide up the spectrum: one wide band, a few wide bands, or many narrow ones. That choice sets what you can measure. The same scene as panchromatic, multispectral, and hyperspectral imagery (the last reduced to three components for display), left to right. ### Panchromatic A panchromatic sensor records one wide band spanning most of the visible and near-infrared. A single wide band collects much more light than a narrow one, so each pixel can be made smaller while still getting a clean signal. That's why the panchromatic band is usually the sharpest image a satellite produces. What it gives up is colour: a panchromatic image records only how bright each pixel is. A panchromatic band: fine detail, brightness only. Panchromatic imagery is often combined with lower-resolution colour bands from the same satellite in a process called pansharpening, which uses the panchromatic detail to sharpen the colour image. It works because the panchromatic band covers the same wavelengths as the colour bands it's sharpening. :::warning Only pansharpen bands that fall inside the panchromatic band's wavelength range. Outside it, the detail you're transferring doesn't correspond to what the band measured, and the sharpened values are unreliable. ::: ### Multispectral A multispectral sensor records a handful of wide bands, typically blue, green, red, and near-infrared. Some add red edge bands, between red and near-infrared, and SWIR bands. Sentinel-2 has 13 bands and Landsat 8 and 9 have 11. Wide bands can overlap, so light at one wavelength registers in two bands at once. Figure 1 shows this for a typical four-band sensor, where the blue and green bands overlap considerably. Figure 1: Spectral response of a typical four-band multispectral satellite. Where two curves overlap, both bands record the same light. `[1]` Figure 2 shows reflectance spectra for soil, vegetation, and water, with typical multispectral bands as grey bars. The bands leave most of the spectrum unsampled, but they still separate these three materials easily, because vegetation reflects far more near-infrared light than soil or water. Figure 2: Soil, vegetation, and water spectra, with typical multispectral bands in grey. Three very different shapes, easy to tell apart even with a few wide bands. `[2]` That works when the materials are very different. It stops working when they're similar. ### Hyperspectral A hyperspectral sensor records many narrow, adjacent bands. Instead of one number averaged over a 100 nm wide band, you get a sample every 10 to 20 nm or so, fine enough to resolve individual absorption features: the chlorophyll absorption near 670 nm, the steep rise of the red edge, or the narrow dips that identify a particular mineral. Figure 3: Hyperspectral band spacing across the VNIR. Each band is narrow, and together they sample the spectrum continuously. Figure 4 shows why that matters. These vegetation types have spectra with the same overall shape, and they differ in small, specific places. Narrow bands resolve those differences. Wide multispectral bands average across them, which makes the types much harder to tell apart. The same applies to soils, minerals, and water constituents. Figure 4: Reflectance spectra of several vegetation types. The differences between them are narrower than a multispectral band. `[3]` Wyvern's Dragonette satellites image the VNIR from 445 to 870 nm, in 23 or 31 bands depending on the satellite. That range covers plant pigments, the red edge, water colour, and iron-bearing minerals. The SWIR, beyond 870 nm, holds features our sensors don't reach, such as leaf and soil water, clays, and carbonates; see [additional datasets](../more/additional_datasets) for missions that cover it. More bands also means bigger files. A hyperspectral scene is many times the size of a multispectral one over the same area, which affects storage and processing time. Our tutorials show how to read only the part of a scene you need. ## Where to go next - [Hyperspectral signatures](hyperspectral_signatures): how a pixel's spectrum identifies what's in it, and what happens when a pixel contains more than one material. - [Plant health indices](../datasets_and_use_cases/plant_health_indices): the most common use of narrow VNIR bands. - [Loading data in Python](/docs/tutorials/python/loading_data): open a real Dragonette scene and look at its bands. ## References `[1]` Schill, Steven & McNulty, Valerie Pietsch & Pollock, F. & Lüthje, Fritjof & Li, Jiwei & Knapp, David & Kington, Joe & McDonald, Trevor & Raber, George & Escovar-Fadul, Ximena & Asner, Gregory. (2021). Regional High-Resolution Benthic Habitat Data from Planet Dove Imagery for Conservation Decision-Making and Marine Planning. Remote Sensing. 13. 4215. 10.3390/rs13214215. `[2]` Pham, Binh. (2018). Satellite remote sensing of the variability of the continental hydrology cycle in the lower Mekong basin over the last two decades. 10.13140/RG.2.2.27136.07686. `[3]` Govender, Megandhren & Chetty, Kershani & Bulcock, Hartley. (2007). A review of hyperspectral remote sensing and its application in vegetation and water resource studies. Water S.A. 33. 10.4314/wsa.v33i2.49049. `[4]` Ronan, Gringer. EM spectrum, revised diagram with re-aligned spectrum. `commons.wikimedia.org/wiki/File:EM_spectrumrevised.png` --- [Source: https://knowledge.wyvern.space/docs/documentation/more/additional_datasets] # Additional datasets Dragonette images the VNIR from 445 to 870 nm at 5.3 m. When you need wavelengths beyond that, a longer archive, or a different scale, these are the datasets to reach for. All three cover the SWIR as well as the VNIR, and all are free to use, though some ask you to register. ## AVIRIS The Airborne Visible / Infrared Imaging Spectrometer (AVIRIS) is NASA JPL's airborne imaging spectrometer. It flies on aircraft rather than satellites, so its pixel size depends on flight altitude, and it covers the full VNIR and SWIR range. It has been flying since 1986 and has built one of the longest-running hyperspectral archives available, used for mineral mapping, vegetation and water studies, and disaster response. There have been three generations of the instrument: AVIRIS-Classic, AVIRIS-NG (Next Generation), and AVIRIS-3, which is the one currently flying. Each has its own data portal: [AVIRIS Data Portal](https://aviris.jpl.nasa.gov/dataportal/): imagery from the original AVIRIS-Classic sensor. This portal holds the most data, since it covers roughly 15 years of flights, though coverage is mostly Canada and the USA. [AVIRIS-NG Data Portal](https://avirisng.jpl.nasa.gov/dataportal/): imagery from the AVIRIS-NG sensor. AVIRIS-NG overflights include Europe, India, and the Arctic. For AVIRIS-3, and for the current list of campaigns and portals, start from the [AVIRIS project site](https://aviris.jpl.nasa.gov/). ## EnMAP The [Environmental Mapping and Analysis Program](https://www.enmap.org/) is a German hyperspectral satellite. It delivers imagery at 30 m GSD in 230 bands from 420 to 2,450 nm, covering the VNIR and SWIR. It launched in April 2022 and started routine operations in November 2022. You need to register to access its data. Its 30 m pixels are coarser than Dragonette's 5.3 m, but it covers the SWIR, which Wyvern's VNIR sensors do not. ## NEON The [National Ecological Observatory Network](https://www.neonscience.org/) is run by the US National Science Foundation to collect ecological data across the United States. Among much else, its airborne program collects [high-resolution hyperspectral imagery](https://www.neonscience.org/data-collection/airborne-remote-sensing) across the VNIR and SWIR over the same field sites year after year, alongside LiDAR and ground measurements, which makes it useful for validating what you see from orbit. --- [Source: https://knowledge.wyvern.space/docs/tutorials/python/visualizing_spectra] # Visualizing Spectra & Plotting One of the outputs from this tutorial, a pixels spectral curve! Try it yourself Runs free in your browser with no install or setup, and pulls the imagery straight from our Open Data Program. Or view the notebook source on GitHub Open in Colab This page details how to produce 2-D scatter plots and spectral curve plots with Wyvern data in Python. These plots can be used to understand the various materials, objects, and species of vegetation that are present in an image. ## Getting started The notebook is self-contained and installs its own dependencies, so there is nothing to set up in advance. - **Google Colab (easiest):** open the notebook using the "Open in Colab" badge above, then run the cells top to bottom. - **Local Jupyter:** open the notebook in JupyterLab or VS Code (Python 3.10+) and run it. Either way, the first cell installs the packages that are not already present: ```python %pip install -q numpy matplotlib requests rasterio ``` ## Image we're using We're using the image also used in [Loading Data in Python](loading_data). Check out the renderings of the image there, or on the [Wyvern Open Data Catalog](https://opendata.wyvern.space/#/year/2024/wyvern_dragonette-001_20240930T070744_08fd7f5a_l2a/wyvern_dragonette-001_20240930T070744_08fd7f5a_l2a.json)! ### Spectral Plots These plots show an individual pixel's surface reflectance across the electromagnetic spectrum that Wyvern's imagery senses. The spectral curve can be used to visually identify the pixel, and can be used by machine learning (and other methods!) to classify the spectral curve as a particular material/species of vegetation. We will dive into classification methods in another tutorial. ### Scatter Plots 2-dimensional scatter plots showcase the distribution of pixels across different bands. This can be used to group spectrally similar pixels, or identify where a particular pixel lies on the plot. ## Code First things first, let's import all of our packages. ```python import json import numpy as np import rasterio import matplotlib.pyplot as plt from matplotlib.patches import Rectangle ``` ## Download the Tutorial Files Let's download the GeoTIFF and its STAC JSON metadata! Run the cell below to download these files. Once the download is complete, they'll appear in the same folder as your notebook! ```python import requests from pathlib import Path IMAGE_URL = "https://wyvern-data.com/wyvern_dragonette-001_20240930T070744_08fd7f5a_l2a/wyvern_dragonette-001_20240930T070744_08fd7f5a_l2a.tiff" JSON_URL = "https://wyvern-odp.com/wyvern_dragonette-001_20240930T070744_08fd7f5a_l2a/wyvern_dragonette-001_20240930T070744_08fd7f5a_l2a.json" IMAGE_PATH = Path("wyvern_dragonette-001_20240930T070744_08fd7f5a_l2a.tiff") JSON_PATH = Path("wyvern_dragonette-001_20240930T070744_08fd7f5a_l2a.json") downloads = [ (IMAGE_URL, IMAGE_PATH), (JSON_URL, JSON_PATH), ] for url, path in downloads: if path.exists(): print(f"{path.name} File already here! Skipping download.") else: print(f"Downloading {path.name}...") response = requests.get(url, headers={"User-Agent": "Mozilla/5.0"}) response.raise_for_status() path.write_bytes(response.content) print("Downloads complete!") ``` ### Loading the Data We'll start by reading the scale factor from the metadata JSON file we downloaded earlier. Wyvern's L2A products represent surface reflectance, but the values are stored as scaled 16-bit integers to reduce file size. ```python # Open the JSON metadata file and extract the scale factor with open(JSON_PATH) as json_file: metadata = json.load(json_file) # The scale factor converts stored uint16 reflectance values into their floating-point representation scale_factor = metadata["assets"]["Cloud optimized GeoTiff"]["raster:bands"][0]["scale"] # if you don't have the JSON file, you can hardcode the scale factor: # scale_factor = 0.0001 print(f"Scale factor: {scale_factor}") ``` Now let's open the GeoTIFF using [rasterio](https://rasterio.readthedocs.io/en/stable/api/index.html). We'll replace any nodata pixels with `np.nan` so they are excluded from our calculations, then apply the scale factor to recover the floating-point reflectance values. ```python with rasterio.open(IMAGE_PATH) as image_file: # Convert to float so we can use np.nan and decimal values image_arr = image_file.read().astype(np.float32) nodata = image_file.nodata band_descriptions = image_file.descriptions # Replace nodata pixels with NaN so they are excluded from calculations if nodata is not None: image_arr[image_arr == nodata] = np.nan # Apply the scale factor to recover floating-point reflectance values image_arr *= scale_factor print(f"Image shape: {image_arr.shape}") print(f"Image bands: {', '.join(band_descriptions)}") ``` ## Building the Wavelength Map Just like in [Loading Data in Python](loading_data), we'll build a map from wavelength to band index using the `eo:bands` metadata, and a small helper to select bands by wavelength. ```python with open(JSON_PATH) as json_file: metadata = json.load(json_file) wavelength_map = {} for i, band in enumerate(metadata["assets"]["Cloud optimized GeoTiff"]["eo:bands"]): wavelength = round(band["center_wavelength"] * 1000) # convert from μm to nm wavelength_map[wavelength] = i print(wavelength_map) ``` ### Selecting Bands by Wavelength ```python def get_band(wavelength_nm: int, return_index: bool = False): """ Returns the image band closest to the requested wavelength. Args: wavelength_nm: Target wavelength in nanometres. return_index: If True, also return the band's index alongside the array. Returns: A 2D array containing reflectance values for the closest available wavelength band, or a (array, band_index) tuple if return_index is True. """ available_wavelengths = list(wavelength_map.keys()) closest_wavelength = min(available_wavelengths, key=lambda w: abs(w - wavelength_nm)) band_index = wavelength_map[closest_wavelength] print(f"Requested {wavelength_nm}nm -- using {closest_wavelength}nm") if return_index: return image_arr[band_index], band_index return image_arr[band_index] ``` # Spectral Curve Plot Now, we're ready to plot our spectral curve. We'll select an arbitrary row/column coordinate in the image to select and plot. Feel free to change these values and explore different pixels in the image! We'll use matplotlib's [subplots](https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.subplots.html) function to easily plot two things in one plot. This will allow us to show the selected pixel's spectra, along with where that pixel lies in the image! This will help us to provide more context on the spectral curve we are looking at. ```python # To plot a spectral curve, we first need to pick a pixel in the image to plot! PIXEL_ROW, PIXEL_COL = 1450, 1800 # Point on image to plot, feel free to change! pixel_spectra = image_arr[:, PIXEL_ROW, PIXEL_COL] # Grabbing pixel from image array # Order our wavelengths by band index so they line up with pixel_spectra for the line plot idx_to_wavelength = {i: wavelength for wavelength, i in wavelength_map.items()} band_wavelengths = [idx_to_wavelength[i] for i in range(image_arr.shape[0])] # Subplots allows us to plot two things at once! Neat! fig, axs = plt.subplots(1, 2, figsize=(13, 5), width_ratios=[1.5, 1]) # First, plotting the spectra as a line chart axs[0].plot(band_wavelengths, pixel_spectra, marker="o", color="#ff4848") axs[0].set_title("Pixel Spectral Curve") axs[0].set_xlabel("Band Wavelength (nm)") axs[0].set_ylabel("Surface Reflectance") # Secondly, plotting a buffered area around our pixel to visualize what we're looking at! rgb_crop = np.dstack([ get_band(650)[PIXEL_ROW - 200 : PIXEL_ROW + 200, PIXEL_COL - 200 : PIXEL_COL + 200], get_band(550)[PIXEL_ROW - 200 : PIXEL_ROW + 200, PIXEL_COL - 200 : PIXEL_COL + 200], get_band(503)[PIXEL_ROW - 200 : PIXEL_ROW + 200, PIXEL_COL - 200 : PIXEL_COL + 200], ]).astype(np.float32) brightness_factor = 3 # Bring up the brightness a touch so we have a better visual rgb_crop_bright = rgb_crop * brightness_factor rgb_crop_bright = np.clip(rgb_crop_bright, 0, 1) # keep values within 0-1 axs[1].imshow(rgb_crop_bright) # We can also add a red box around the pixel we are looking at in the line chart! pixel_rect = Rectangle( (195, 195), 10, 10, linewidth=2, edgecolor="#ff4848", facecolor="none" ) axs[1].add_patch(pixel_rect) axs[1].set_title("Pixel Location in Image (RGB)") plt.tight_layout() plt.show() ``` Spectral curve of selected pixel. Red box denotes where the pixel resides. ## Spectral Curve Discussion This pixel shows there is vegetation present. The reflectance stays low and fairly flat through the visible range, since chlorophyll is absorbing most of that light, then jumps sharply around 680–730nm into the near-infrared. That jump is called the "red edge," and this makes healthy vegetation easy to pick out from soil or other surfaces! Because this is L2A surface reflectance data, atmospheric effects have already been largely corrected for, so the curve above is an estimate of surface reflectance after atmospheric correction, much closer to the material's actual reflectance at this pixel than raw radiance would be. Keep in mind that comparing it against reference spectral libraries still requires accounting for differences in available bands and bandwidths between Wyvern's sensor and the library measurements, which we cover in more detail in the [Detecting Rare Earth Elements](detecting_rare_earth_elements) tutorial. # Scatter Plot/Heatmap Heatmap/scatter plots of pixel distributions are useful for discerning different groups of pixels within an image. For example, one group of pixels may be vegetation, while another group of pixels might be dirt or asphalt. ```python # Now, we're ready to plot some band combinations! # We'll start out w/ Bands 649nm and 764nm band_x, X_IDX = get_band(649, return_index=True) band_y, Y_IDX = get_band(764, return_index=True) # Ravel turns our 2D X-Y image array into a single dimensional array band_x = band_x.ravel() band_y = band_y.ravel() # Here we're filtering out our NaNs from the array. We use a single joint mask so that if a # pixel is NaN in either band, it's dropped from both, keeping band_x and band_y aligned pixel-for-pixel. valid = ~np.isnan(band_x) & ~np.isnan(band_y) band_x = band_x[valid] band_y = band_y[valid] # This generates a two-dimensional histogram, or summary of our data statistics. # It allows us to plot a "heatmap" of where points lie between two bands. # We use percentiles to set the plot visualization range. This just # keeps a few very bright or very dark pixels from stretching the axes x_min, x_max = np.percentile(band_x, [0.01, 99.9]) y_min, y_max = np.percentile(band_y, [0.01, 99.9]) heatmap, xedges, yedges = np.histogram2d( band_x, band_y, bins=200, range=[[x_min, x_max], [y_min, y_max]] ) extent = [xedges[0], xedges[-1], yedges[0], yedges[-1]] # We then can plot our heatmap! fig, ax = plt.subplots(1, 1) ax.imshow(heatmap.T, extent=extent, origin="lower", cmap="viridis") plt.xlabel(band_descriptions[X_IDX]) plt.ylabel(band_descriptions[Y_IDX]) # Since we highlighted where our selected pixel was on the image in the above spectral plot, let's # also highlight where the pixel lies in our scatter plot! box_size = 0.015 pixel_rect = Rectangle( (pixel_spectra[X_IDX] - box_size / 2, pixel_spectra[Y_IDX] - box_size / 2), box_size, box_size, linewidth=2, edgecolor="#ff4848", facecolor="none" ) ax.add_patch(pixel_rect) plt.title("2 Dimensional Band Scatter Plot") plt.tight_layout() plt.show() ``` 649nm & 764nm Scatter Plot. Red box denotes pixel we looked at in the spectra plots. ```python # Let's now do our first and last bands! band_x, X_IDX = get_band(503, return_index=True) band_y, Y_IDX = get_band(799, return_index=True) # Ravel turns our 2D X-Y image array into a single dimensional array band_x = band_x.ravel() band_y = band_y.ravel() # Here we're filtering out our NaNs from the array. We use a single joint mask so that if a # pixel is NaN in either band, it's dropped from both, keeping band_x and band_y aligned pixel-for-pixel. valid = ~np.isnan(band_x) & ~np.isnan(band_y) band_x = band_x[valid] band_y = band_y[valid] # This generates a two-dimensional histogram, or summary of our data statistics. # It allows us to plot a "heatmap" of where points lie between two bands. # We use percentiles to set the plot visualization range. This just # keeps a few very bright or very dark pixels from stretching the axes x_min, x_max = np.percentile(band_x, [0.01, 99.9]) y_min, y_max = np.percentile(band_y, [0.01, 99.9]) heatmap, xedges, yedges = np.histogram2d( band_x, band_y, bins=200, range=[[x_min, x_max], [y_min, y_max]] ) extent = [xedges[0], xedges[-1], yedges[0], yedges[-1]] # We then can plot our heatmap! fig, ax = plt.subplots(1, 1) ax.imshow(heatmap.T, extent=extent, origin="lower", cmap="viridis") plt.xlabel(band_descriptions[X_IDX]) plt.ylabel(band_descriptions[Y_IDX]) # Since we highlighted where our selected pixel was on the image in the above spectral plot, let's # also highlight where the pixel lies in our scatter plot! box_size = 0.015 pixel_rect = Rectangle( (pixel_spectra[X_IDX] - box_size / 2, pixel_spectra[Y_IDX] - box_size / 2), box_size, box_size, linewidth=2, edgecolor="#ff4848", facecolor="none" ) ax.add_patch(pixel_rect) plt.title("2 Dimensional Band Scatter Plot") plt.tight_layout() plt.show() ``` 503nm & 799nm Scatter Plot. Red box denotes pixel we looked at in the spectra plots. As-generated these plots do not tell us much on their own. In both plots we can see a tight, bright cluster of pixels, this is the dirt/sand around the irrigated fields. The Y-axis in both plots represents a near-infrared band, while the X-axis represents a visible band. Vegetation generally has high near-infrared reflectance and lower visible reflectance, which produces the fainter diagonal spread of pixels extending away from that cluster toward higher Y (and lower X) values. We can see some distinct groups within that spread (feel free to zoom the plot in to this area!), these could correspond to different crop types, growth stages, or vegetation conditions. How cool! Our vegetation pixel lands in the top left corner of both plots. Plants absorb most visible light for photosynthesis, keeping their visible band reflectance low, while reflecting near-infrared light, pushing their NIR band reflectance high. That combination of low visible and high NIR reflectance is what puts vegetation in that top left region. We could extend these plots by highlighting areas of high-NDVI, or labeling areas with a shapefile and highlighting these areas to identify pixels that have a matching/similar spectra. We'll leave extending this notebook up to you, dear reader :) --- [Source: https://knowledge.wyvern.space/docs/tutorials/qgis/pca] # Creating PCA images in QGIS Image of Saudi Arabia agriculture: Natural colour on left, 3-band PCA on right In this tutorial, we're going to use QGIS to apply Principal Component Analysis (PCA) to a Wyvern hyperspectral image. PCA is a form of [unsupervised learning](https://en.wikipedia.org/wiki/Unsupervised_learning). When applied to hyperspectral imagery, it creates a new set of bands (or components) which are: - unrelated to each other (orthogonal) - and ordered by the amount of image variation they explain Because the bands are unrelated to each other, each band contains unique information about the original image. And because the bands are ordered by the image variation they explain, the first few bands will contain the most information about the image; you can focus your attention on those. The end result is an image that can be used to reveal hidden correlations and similarities within the original image, and highlight areas for further exploration. ## Before you start... Before starting this tutorial, please ensure you have done the following: - QGIS has been downloaded installed on your machine. Visit [qgis.org](https://www.qgis.org/) to download QGIS. - You have [installed the Orfeo Toolbox plugin into QGIS](https://www.orfeo-toolbox.org/CookBook/QGISInterface.html). - You have a Wyvern image you can load into QGIS. If you haven't downloaded one already, we'll be [using this image](https://opendata.wyvern.space/#/year/2024/wyvern_dragonette-001_20240930T070744_08fd7f5a/wyvern_dragonette-001_20240930T070744_08fd7f5a.json) for the next steps. :::info **Note:** Apple silicon devices may have trouble installing QGIS. For best results, use a Windows or Linux (x86) machine, or follow the instructions on the QGIS website to install the OSX version. ::: ## Generating PCA images in QGIS :::info **Note:** OTB and `NODATA` values OTB does not currently have a simple way to handle `NODATA` values. Because of this, we're going to crop our raster to exclude `NODATA` values, and perform PCA on that. More info can be found [in this OTB issue](https://gitlab.orfeo-toolbox.org/orfeotoolbox/otb/-/issues/2441). ::: ### Cropping the image Ensure the image is selected in the `Layers` pane. Click on the `Raster` menu, then `Extraction`, then `Clip raster by Extent`. In the dialog box, click on the `Clipping extent` menu and select `Draw on map canvas`. Clip raster by Extent dialog Draw on the map, being careful to stay within the borders of the displayed layer, then click `Run`. You should now see a new layer named `Clipped (extent)` with all 23 bands: Clipped (extent) layer ### Generate the PCA Click on the Processing menu, then Toolbox. In the "Search..." bar, type `DimensionalityReduction`. You should see it listed under `OTB` -> `Image Filtering`. QGIS Dimensionality Reduction Menu Click on that, and you should see a new dialog window pop up: QGIS Dimensionality Reduction Dialog - For `Input Image`, ensure the `Clipped (extent)` layer is selected. - For `Output image`, click on the 3 dots and select the directory where you're keeping the original image. Otherwise, it will use a temporary file, and you may run out temp space depending on your computer. Name the file `clipped_pca.tif` for now. - For `Inverse Output image`, `Output files containing eigenvalues` and `Transformation matrix`, select `Skip output`. Click `Run`. Within a minute or two, the processing should be done and you should see a new layer in the `Layers` window. Now we can start to examine the data! ## Different ways of viewing the PCA images ### RGB image By default, QGIS will select the first three bands (ie, the first three components, or eigenrasters) as the Red, Green and Blue bands. We've set the `Min/ Max Value Settings` to `Mean +/- standard deviation` to make the layer distinct, but you can experiment with these to see what works best for you. PCA Multiband Colour You can also experiment with selecting different bands (ie, different eigenvectors) to see what this brings out in the image. Note that as you choose higher and higher bands, the effects will decrease -- you're effectively looking at lesser and lesser significant sources of information within the image. ### Singleband pseudocolour As with the RGB preview, picking the correct single band pseudocolour for your application will require some exploration. Here, we're showing Band 1 (ie, the most significant eigenraster) in the `Turbo` colour ramp: PCA Single Band Pseudocolour We can see a clear differentiation between not just the irrigated plots (blue), non-irrigated plots (yellow), and the empty plots (green), but the buildings showing up as smaller rectangular areas (red) within the plots at the top. ## What does PCA show us? Now that we've found helpful ways to display the PCA bands, let's dig into what PCA might tell us. We'll focus on a few areas. ### Striking similarities between different sections of the image Here we have a zoomed-in area showing fields, roads and paths: Natural colour PCA - Bands 1, 2, 3 Roads and fields: Natural colour vs PCA Bands 1, 2, 3 We'll highlight a few areas of similarity here from the PCA image: - The yellow areas in the PCA image consist of field sections and borders, roads, and irregular paths. - The purple border of the field in the top left quadrant matches field sections in the bottom right. - The greens in the three fields in the top left quadrant are quite close, despite being readily distinguishable in natural colour. If we compare these parts of the PCA image to the natural colour image, these correspondences are either not apparent or not nearly as striking. They suggest a relationship between disparate areas that are not immediately visible in the original natural colour image. ### Differences within similarly-coloured areas are much more apparent Here, we take a closer look at four adjacent fields. Pay attention to the one in the bottom left: Natural colour PCA - Bands 1, 2, 3 Roads and fields: Natural colour vs PCA Bands 1, 2, 3 PCA has brought out differences within that field, particularly the structure in the top half, that are not nearly as obvious in the original natural colour image. ### Different PCA bands contain different information Finally, we take a look at another set of fields -- but this time, we are changing the bands we use. Below on the left, you can see the image composed of PCA Bands 1, 2 and 3; on the right, it's composed of PCA Bands 2, 3, and 4: PCA - Bands 1, 2, 3 PCA - Bands 2, 3, 4 PCA: Bands 1, 2, 3 vs Bands 2, 3, 4 Here are a few things the change in bands bring out: - The structure of the borders between field sections, and the outer borders of the fields - Additional detail in the centre of each circular field - We *lose* some features in the rectangular field on the right but highlight others This highlights the value of exploring different bands to see what each one reveals. :::info **Note** PCA bands are ordered by how much information they contain about the original image. A common technique for simplifying PCA results is choosing to discard some of the lower bands, as they're "less important" than the higher ones. Choosing where to draw that line is often subjective, and outside the scope of this tutorial -- but we'll point you to [scree plots](https://en.wikipedia.org/wiki/Scree_plot) and the [elbow method](https://en.wikipedia.org/wiki/Elbow_method_(clustering)) to get started. ::: ## Next steps We've focused here on one image (Saudia Arabia agriculture), and on the first three components produced by PCA -- but that leaves much more to explore. Here are some questions to investigate: - What does the RGB image look like if you select a different set of bands (components) from the PCA layer? What differences are highlighted? - What happens as you select higher and higher bands from the PCA layer? What does this tell you about the image? How would you decide which bands to ignore? - Try looking at pseudocolour images of individual bands. How is this useful compared to an RGB image? - What happens if you pick another image from our [open data program](https://opendata.wyvern.space)? What does PCA highlight if you switch from agriculture to [urban](https://opendata.wyvern.space/#/wyvern_dragonette-002_20250503T051751_81db60aa/wyvern_dragonette-002_20250503T051751_81db60aa.json), [forest](https://opendata.wyvern.space/#/wyvern_dragonette-002_20250728T231028_6da86460/wyvern_dragonette-002_20250728T231028_6da86460.json), [mining](https://opendata.wyvern.space/#/wyvern_dragonette-001_20240811T083914_08c61457/wyvern_dragonette-001_20240811T083914_08c61457.json), or [maritime](https://opendata.wyvern.space/#/year/2024/wyvern_dragonette-001_20240628T062939_5fce57a3/wyvern_dragonette-001_20240628T062939_5fce57a3.json) scenes? ## Troubleshooting ### The first band is all one colour! This can happen if you run PCA on the entire original image, rather than selecting an extent and running PCA on that. The reason is that our image contains `NODATA` values surrounding the image; when PCA sees that, it decides that the first and most important way the data varies is by whether it's `NODATA` or not. ### I get an error saying there was a problem processing the image! Symptoms: a message that says `The following layers were not correctly generated`. Check your tmp directory to see if it's filled up -- this may be left over from previous attempts to generate the PCA image. You should be able to remove the files *within* the directory without problems, and try again. You can also try setting a different tmpdir for this process. Click on Advanced, then Parameters, and set `Temporary Folder` as needed. --- [Source: https://knowledge.wyvern.space/docs/documentation/getting_started/hyperspectral_image_formats] # Hyperspectral image formats A hyperspectral image is a raster like any other, just with far more bands. This page explains how to think about that structure and covers the three file formats you're most likely to meet: GeoTIFF (which Wyvern delivers), ENVI, and HDF5. The number of bands is the main thing that sets hyperspectral data apart, and it's also why the files are large: | Satellite | Type | Bands | Common format | | --------- | ---- | ----- | ------------- | | Sentinel-2 | Multispectral | 13 | JPEG2000, GeoTIFF | | Landsat 8 and 9 | Multispectral | 11 | GeoTIFF | | Sentinel-1 | SAR | 1–2 polarizations | GeoTIFF | | **Wyvern Dragonette**| **Hyperspectral** | **23 or 31** | **Cloud-Optimized GeoTIFF** | ## The data cube A hyperspectral scene is easiest to picture as a cube: two spatial axes, x and y, and one spectral axis, the bands. Slice the cube at a single band and you get a greyscale image. Take three slices and assign them to red, green, and blue, and you get a colour composite. Take one pixel and read straight down through every band, and you get that pixel's spectrum, which is where most hyperspectral analysis starts. A hyperspectral data cube. The top face is the image; depth is wavelength. `[1]` Once it's loaded, the file format stops mattering: every raster becomes an array of numbers with a location attached. If you're not sure whether a tool can read a file, try GDAL. It reads almost everything. ## Common formats These are the three formats you'll meet most often with hyperspectral data. GDAL's [raster driver list](https://gdal.org/drivers/raster/index.html) covers the rest. ### GeoTIFF `.tiff` GeoTIFF is a TIFF image with embedded georeferencing, so a single file holds both the pixels and where they sit on the Earth. It was developed by Dr. Niles Ritter at NASA JPL and is now an Open Geospatial Consortium standard. It has variants for different settings, including BigTIFF for very large files and Cloud-Optimized GeoTIFF. Wyvern delivers imagery as Cloud-Optimized GeoTIFF (COG). A COG is a regular GeoTIFF with its data arranged in tiles and overviews, so software can read just the part it needs over HTTP instead of downloading the whole file. Each band's centre wavelength and FWHM are stored as tags in the file, so you get the spectral metadata along with the pixels. [Wikipedia article on the GeoTIFF standard](https://en.wikipedia.org/wiki/GeoTIFF) [Open Geospatial Consortium GeoTIFF standard](https://www.ogc.org/standards/geotiff) [GDAL GeoTIFF raster driver](https://gdal.org/drivers/raster/gtiff.html) ### ENVI `.img` and `.hdr` The ENVI format is the native format of NV5's ENVI software, and it's widely used for hyperspectral data. It comes in two parts: a binary image file and a plain-text header, usually with a `.hdr` extension. The header is where ENVI stores band names, wavelengths, and FWHM, and because it's plain text you can open it in any editor to check what a file contains. AVIRIS data is distributed in this format. GDAL (along with rasterio and any other tool based on GDAL) can read and write it. [NV5 ENVI documentation](https://www.nv5geospatialsoftware.com/docs/enviimagefiles.html) [GDAL ENVI raster driver](https://gdal.org/drivers/raster/envi.html) ### HDF5 `.h5` HDF5 is a scientific data format maintained by the non-profit HDF Group. It's organized like a file system: named groups contain datasets or further groups, and metadata is attached to each. NetCDF4 is built on HDF5 and is worth knowing about because it's what NASA's EMIT mission uses, and because the `xarray` library in Python reads it directly, keeping band and coordinate labels attached to the array. PRISMA uses an HDF5 variant, HDF-EOS5. [Wikipedia article on the HDF standard](https://en.wikipedia.org/wiki/Hierarchical_Data_Format) [GDAL HDF5 raster driver](https://gdal.org/drivers/raster/hdf5.html) ## More resources [GIS Geography's list of GIS formats](https://gisgeography.com/gis-formats/) [GDAL raster driver list](https://gdal.org/drivers/raster/index.html) ## References `[1]` Pizzolante, Raffaele & Carpentieri, Bruno. (2012). Visualization, Band Ordering and Compression of Hyperspectral Images. Algorithms. 5. 10.3390/a5010076. --- [Source: https://knowledge.wyvern.space/docs/tutorials/python/index_calculation] # Index Calculations with Python Calculating indices like NDVI or NDWI is an extremely common task in geospatial analysis. Many tools like QGIS provide a user-friendly interface for calculating indices; with Python we have to do more of the hard work ourselves. Using a more programmatic interface for these analyses lets us start to build data pipelines and scale our data solutions. Red-Edge NDVI (RENDVI) that we'll render in this tutorial from Wyvern L2A hyperspectral imagery over irrigation fields in southern Alberta. Index calculations, in their simplest form, is just math applied to arrays of numbers (raster images!). There are a number of packages available in Python for working with arrays. We'll cover the most popular (Numpy) here. This tutorial will use this image over southern Alberta, available for download at our [Open Data Program](https://opendata.wyvern.space/#/year/2025/wyvern_dragonette-004_20251005T172442_932fa466_l2a/wyvern_dragonette-004_20251005T172442_932fa466_l2a.json). --- ## High Level Workflow 1. Load in our raster data using [Rasterio](https://rasterio.readthedocs.io/en/stable/). 2. Create a dictionary the band wavelengths and their corresponding index in the raster. 3. Calculate our indices using [Numpy](https://numpy.org/). Try it yourself Runs free in your browser with no install or setup. Or view the notebook source on GitHub Open in Colab --- ## Getting started The companion notebook is self-contained and installs its own dependencies, so there is nothing to set up in advance. - **Google Colab (easiest):** open the notebook using the "Open in Colab" badge at the top of it, then run the cells top to bottom. - **Local Jupyter:** open the notebook in JupyterLab or VS Code (Python 3.10+) and run it. Either way, the first cell installs the packages that are not already present: First let's import all the packages we'll need. You may need to install them using `pip install -q numpy matplotlib requests raster`. ```python import rasterio import numpy as np import json import matplotlib.pyplot as plt ``` ## Download the Tutorial Files Let's download the files we'll be working with for this tutorial! We'll be using a hyperspectral GeoTIFF image over southern Alberta and a STAC JSON metadata file. The GeoTIFF contains the image data, and the STAC JSON gives us useful metadata about the scene. Run the cell below to download these files. Once the download is complete, they'll appear in the same folder as your notebook! ```python import requests from pathlib import Path IMAGE_URL = "https://wyvern-data.com/wyvern_dragonette-004_20251005T172442_932fa466_l2a/wyvern_dragonette-004_20251005T172442_932fa466_l2a.tiff" JSON_URL = "https://wyvern-odp.com/wyvern_dragonette-004_20251005T172442_932fa466_l2a/wyvern_dragonette-004_20251005T172442_932fa466_l2a.json" IMAGE_PATH = Path("wyvern_dragonette-004_20251005T172442_932fa466_l2a.tiff") JSON_PATH = Path("wyvern_dragonette-004_20251005T172442_932fa466_l2a.json") downloads = [ (IMAGE_URL, IMAGE_PATH), (JSON_URL, JSON_PATH), ] for url, path in downloads: if path.exists(): print(f"{path.name} File already here! Skipping download.") else: print(f"Downloading {path.name}...") response = requests.get(url, headers={"User-Agent": "Mozilla/5.0"}) response.raise_for_status() path.write_bytes(response.content) print("Downloads complete!") ``` ### Loading the Data We'll start by reading the scale factor from the metadata JSON file we downloaded earlier. Wyvern's L2A products represent surface reflectance, but the values are stored as scaled 16-bit integers to reduce file size. We can grab the scale value we need from the metadata file. More documentation about processing levels can be found in our [documentation](https://guide.wyvern.space/data-product-guide/technical-specifications/processing-levels/#level-2a-surface-reflectance). ```python # Open the JSON metadata file and extract the scale factor with open(JSON_PATH) as json_file: metadata_json = json.load(json_file) scale_factor = metadata_json["assets"]["Cloud optimized GeoTiff"]["raster:bands"][0]["scale"] print(f"Scale factor: {scale_factor}") ``` Now let's open the GeoTIFF. We'll replace any nodata pixels (the provided uint16 maximum of `65535`) with `np.nan` so they are excluded from our calculations, then apply the scale factor to recover the floating-point reflectance values. ```python with rasterio.open(IMAGE_PATH) as image_file: # Convert to float so we can use np.nan and decimal values image_array = image_file.read().astype(float) nodata_value = image_file.nodata # Arrange nodata values to be np.nan if nodata_value is not None: image_array[image_array == nodata_value] = np.nan # Scale data based on scale factor in metadata file image_array = image_array * scale_factor ``` Now let's build a map from wavelength to band index using the `eo:bands` metadata, so we can select bands by wavelength rather than hardcoding band numbers. This lets our code be more reusable across Wyvern scenes/sensors with different band configurations. ```python eo_bands = metadata_json["assets"]["Cloud optimized GeoTiff"]["eo:bands"] wavelength_map = {} for i, band in enumerate(eo_bands): wavelength = round(band["center_wavelength"] * 1000) # convert from μm to nm wavelength_map[wavelength] = i print(wavelength_map) ``` ### Selecting Bands by Wavelength Now we'll write a small helper function that finds the nearest available band to any wavelength we ask for. This means if we request 550nm and the closest available band is actually 549nm, it'll find it automatically and tell us which band it used. ```python def get_band(wavelength_nm: int) -> np.ndarray: """ Returns the image band closest to the requested wavelength. Args: wavelength_nm: Target wavelength in nanometres. Returns: A 2D array containing reflectance values for the closest available wavelength band. """ available_wavelengths = list(wavelength_map.keys()) closest_wavelength = min(available_wavelengths, key=lambda w: abs(w - wavelength_nm)) print(f"Requested {wavelength_nm}nm -- using {closest_wavelength}nm") return image_array[wavelength_map[closest_wavelength]] ``` Let's take a look at our image and confirm everything was read in correctly. We'll use Matplotlib and plot an RGB rendering of our data. ```python # Select our RGB bands by wavelength, then stack them into a [rows, columns, channels] # array which matplotlib expects rgb_arr = np.dstack([get_band(650), get_band(550), get_band(480)]).astype(np.float32) brightness_factor = 3 # Bring up the brightness a touch so we have a better visual rgb_arr_bright = rgb_arr * brightness_factor rgb_arr_bright = np.clip(rgb_arr_bright, 0, 1) # keep values within 0-1 plt.figure(figsize=(12, 8)) plt.imshow(rgb_arr_bright) plt.axis("off") plt.show() ``` Natural color composite of center pivot irrigation fields in south-eastern Alberta. That all looks great! Now we will move to actually compute the indices of use to us for our image. We will compute RENDVI and MTVI2.\ The wavelengths selected are based off the data available in our [hyperspectral library](../../../hyperspectral_library); feel free to explore for other indices you may be interested in! Note that we're using `get_band()` here, so each of these is the *nearest available* Wyvern band to the nominal index wavelength, not necessarily an exact match. Keep an eye on the printed "Requested Xnm -- using Ynm" output to see how close each one lands. ```python # RENDVI band_750 = get_band(750) band_705 = get_band(705) # MTVI2 band_550 = get_band(550) band_670 = get_band(670) band_800 = get_band(800) # Compute our indices rendvi_array = (band_750 - band_705) / (band_750 + band_705) mtvi2 = 1.5 * ( (1.2 * (band_800 - band_550) - 2.5 * (band_670 - band_550)) / ( np.sqrt( (2 * band_800 + 1) ** 2 - (6 * band_800) + (5 * np.sqrt(band_670)) - 0.5 ) ) ) ``` Now we can plot our RGB image alongside our freshly computed indices and see how they look. ```python # Plot our bright RGB and two indices fig, axes = plt.subplots(1, 3, figsize=(14, 6)) # Plot RGB image im0 = axes[0].imshow(rgb_arr_bright) axes[0].set_title("RGB") axes[0].axis("off") # Plot Red-Edge NDVI im1 = axes[1].imshow(rendvi_array, cmap="RdYlGn") axes[1].set_title("Red-Edge NDVI") axes[1].axis("off") fig.colorbar(im1, ax=axes[1], fraction=0.046, pad=0.04) # Plot MTVI2 im2 = axes[2].imshow(mtvi2, cmap="RdYlGn") axes[2].set_title("MTVI2") axes[2].axis("off") fig.colorbar(im2, ax=axes[2], fraction=0.046, pad=0.04) plt.tight_layout() plt.show() ``` RGB, RENDVI, and MTVI2 indices plotted. Those look great! Hyperspectral imagery's advantage here is its much finer spectral sampling, which lets us select narrower, more precisely-placed bands than a typical multispectral or optical sensor usually offers. Compared to RGB, both hyperspectral indices enhance biophysical signal over visual appearance. RENDVI emphasizes chlorophyll-related spectral shifts in the red-edge region, whereas MTVI2 is designed to stay sensitive to green leaf area and canopy structure while reducing the influence of chlorophyll variation and soil background. The differences in spatial patterns between the two indices demonstrate how band selection and mathematical formulation influence the type of vegetation information extracted from hyperspectral data. Try building these indices on your own dataset and compare the results.\ The more you experiment, the more you’ll uncover. --- ## Other Packages Many Python packages meant for remote sensing data will also include helpful functions for calculating indices. There also exists a wider ecosystem outside of Python for analysis of multispectral & hyperspectral data that is worth diving into! ### EarthPy [EarthPy](https://earthpy.readthedocs.io/en/latest/api/earthpy.spatial.html#earthpy.spatial.normalized_diff) provides a useful `normalized_diff()` function, which can be used to calculate indices such as NDVI ### Xarray-Spatial [Xarray-spatial has built in functions](https://xarray-spatial.org/user_guide/multispectral.html) for calculating a number of useful multispectral indices using Xarray arrays, including NDVI, SAVI, and True Color ### Awesome-Spectral-Indices David Montero Loaiza's [Awesome Spectral Indices](https://github.com/awesome-spectral-indices/awesome-spectral-indices) project contains a number of packages in different languages for calculating multispectral indices. For the Python ecosystem, [Spyndex](https://github.com/awesome-spectral-indices/spyndex) is available for use. The advantage of this system is that a standard database of indices is pulled from, instead of having to form the equations yourself. --- [Source: https://knowledge.wyvern.space/docs/tutorials/qgis/visualizing_spectra] # Visualizing Spectra & Plotting Band 535nm (X) and 799nm (Y) rendered as a 2D density scatter plot in EnMAP-Box ## Before you start... Before starting this tutorial, please ensure you have done the following: - QGIS has been downloaded and installed on your machine. Visit [qgis.org](https://www.qgis.org/) to download QGIS. - You have [installed the EnMAP-Box plugin into QGIS](https://enmap-box.readthedocs.io/en/latest/usr_section/usr_installation.html) - You have loaded a Wyvern image into QGIS and/or EnMAP-Box. If you haven't downloaded one already, we'll be [using this image](https://opendata.wyvern.space/#/year/2024/wyvern_dragonette-001_20240930T070744_08fd7f5a/wyvern_dragonette-001_20240930T070744_08fd7f5a.json) for the next steps. ## Visualizing Spectra in QGIS In the `Attributes` tool bar click the ![](images/identify_features_icon.jpg) `Identify Features` Button, or press the hotkey `Ctrl + Shift + I`. This will open the `Identify Results` pane on the left-hand side of your QGIS window. Click the `View` dropdown and select `Graph` to configure the `Identify Results` pane to display pixel information as a graph. Empty Identify Results pane with the Graph view mode selected Now, you're ready to click on **any part of the image** with the `Identify Features` tool. Once you click on any part of the image, the selected pixel's spectra will show up in the `Identify Results` pane as a line chart! QGIS window with pixel selected & spectra visible ?> Make sure to have the image layer you want to view selected in the `Layers` pane when using the `Identify Features` tool! [More information on the Identify Features tool is available at the QGIS documentation](https://docs.qgis.org/3.34/en/docs/user_manual/introduction/general_tools.html#identify) ## Visualizing Spectra in EnMAP-Box EnMAP-Box has [fantastic documentation on their Spectral Libraries tool](https://enmap-box.readthedocs.io/en/release_3.10/usr_section/usr_manual/spectral_library.html#spectral-libraries) that we will be using. We highly recommend also taking a look at their documentation. To get started, ensure you have opened a Wyvern image in a map view. Then, click on the ![](images/enmap_identify_pixels_button.jpg) `Identify pixel profiles` tool on the top toolbar of the EnMAP-Box application. This will allow us to click on any pixel and get it's spectral curve. It will also automatically open up a new spectral library view for us! When you're ready, click on any pixel to bring up it's spectra in a new spectral library! EnMAP-Box with the selected pixels spectra visible And now we see it's spectra! Well done! If the green on black vibe is not doing it for you, EnMAP-Box has some documentation on how to customize the plot which is [located here](https://enmap-box.readthedocs.io/en/release_3.10/usr_section/usr_manual/spectral_library.html#spectral-profile-sources). If you would like to collect different spectral curves in an image, you can click on the ![](images/enmap_add_spectra_button.jpg) `Add Profile` button to add the selected spectra to the current Spectral Library. You can then view the added spectra (and add a name to the spectra!) by clicking the ![](images/enmap_open_attribute_table.jpg) Open Attribute Table button. This will open the attribute table, which you can then click into to add a name to the spectra you just added! Adding spectra to EnMAP-Box with the Attribute Table open For more information on working with EnMAP-Box's spectral library tool, [take a look at EnMAP-Box's documentation](https://enmap-box.readthedocs.io/en/release_3.10/usr_section/usr_manual/spectral_library.html#spectral-libraries) ## EnMAP-Box Band Scatter Plots Two-dimensional scatter plots allow you to plot two bands against each other, and identify groups of similar pixels. In this example, we will plot bands 535nm and 799nm. To generate 2D band scatter plots, use the top ribbon to select the `Tools -> Scatter Plot` tool. Selecting the Scatter Plot tool This will open the Scatter Plot tool as a separate window. Select your X & Y bands by first selecting the image you have imported into EnMAP-Box, and then selecting the band index for the X and Y values. Your scatter plot will then generate! Rendered scatter plot with Bands 535nm and 799nm Check out the distinct groups in this scatter plot! Move the EnMAP-Box map extent around to see how the scatter plot changes with what pixels are available to plot. You can also adjust the X & Y ranges to cut off outliers. For more information on the Scatter Plot tool, [check out EnMAP-Box's documentation on the tool](https://enmap-box.readthedocs.io/en/latest/usr_section/usr_manual/tools.html#scatter-plot). --- [Source: https://knowledge.wyvern.space/docs/documentation/getting_started/processing_levels] # Processing levels A satellite doesn't measure reflectance. It measures a digital number (DN) for each pixel and band. That number mixes the light reflected by the surface with light absorbed and scattered by the atmosphere, shading from terrain, and the sensor's own response. Turning DNs into something physically meaningful takes several processing steps, and processing levels are the shorthand for how far along that chain a product is. This page explains the standard levels and then the two that Wyvern delivers, L1B and L2A, and how they differ in practice. There are two kinds of correction. *Radiometric* correction converts DNs into physical units: first radiance, then reflectance. *Geometric* correction assembles the image, removes distortions from sensor tilt and terrain, and ties each pixel to a location on the ground. ## Standard processing levels The levels run from 0 to 4. Definitions vary between agencies and sensors, but the framework below is broadly shared. | Level | Description | Wyvern data | | ----------- | ----------- | :-----------: | | Level 0 (L0) | Reconstructed, unprocessed instrument data at full resolution, with communication artifacts removed. Includes the image data, image metadata, and relevant ancillary data. | | | Level 1A (L1A) | Reconstructed, unprocessed instrument data at full resolution, time-referenced and annotated with ancillary information such as camera settings, spacecraft altitude, and orbital ephemeris. | | | Level 1B (L1B) | Image data radiometrically processed to at-sensor, or top-of-atmosphere (TOA), radiance. Includes geolocation and the calibration parameters used to produce it. | :heavy_check_mark: | | Level 2A (L2A) | Image data atmospherically corrected to bottom-of-atmosphere (BOA) surface reflectance. | :heavy_check_mark: | | Level 2B (L2B) | Definitions vary by provider. Commonly, variables derived from L2A reflectance at the same resolution and location. | | | Level 3 (L3) | Data mapped onto a uniform space-time grid, usually averaged over space, time, or both. | | | Level 4 (L4) | Model outputs or analysis results derived from lower-level data, often from several sources. | | ## Wyvern's data products Wyvern delivers two processing levels, and the [Open Data Program](https://opendata.wyvern.space/) publishes both. L2A is the default: the root catalog states that all imagery is L2A surface reflectance unless otherwise noted, and a smaller [L1B collection](https://wyvern-odp.com/l1b/collection.json) holds radiance scenes. Item IDs carry an `_l1b` or `_l2a` suffix, and the STAC item's `processing:level` property states the level directly. The two products are not interchangeable. They differ in data type, scale factor, NoData value, and coordinate reference system, so code written against one will silently produce wrong numbers on the other. | | L1B top-of-atmosphere radiance | L2A surface reflectance | | --- | --- | --- | | Measures | At-sensor radiance | Bottom-of-atmosphere reflectance | | Units | W m-2 sr-1 µm-1 | Unitless reflectance, 0 to 1 | | Data type | `float32` | `uint16` | | Scale factor | None | 0.0001, apply after masking | | NoData | -9999 | 65535 | | Coordinate reference system | Geographic WGS84 ([EPSG:4326](https://epsg.io/4326)) | Per-scene UTM zone, read `proj:epsg` | | Pixel geometry | Non-square in angular degrees | Square metres | | Masks | Data Mask, Pixel Quality Mask | Data Mask, Pixel Quality Mask | Both are delivered as cloud-optimized [GeoTIFF](hyperspectral_image_formats#geotiff-tiff) (COG) alongside [STAC](stac_documentation_and_tools) `.json` metadata, with the metadata you need for further radiometric and geometric processing and for analysis such as spectral indices, classification, anomaly detection, material identification, sub-pixel unmixing, and target detection. :::warning On L2A, mask NoData before you apply the scale factor. Scaling first turns 65535 into 6.5535, a value that looks plausible and corrupts every statistic you compute afterwards. Note also that the GeoTIFF's own `scales` tag is 1.0; the 0.0001 factor comes from the STAC item's `raster:bands`. ::: ### L1B geometry L1B is georeferenced north-up to Geographic WGS84 (EPSG:4326) and radiometrically corrected to at-sensor radiance, using the spacecraft's position and pointing and the solar conditions at the time of acquisition. Because EPSG:4326 is in degrees, pixel size is set to the number of degrees that equals 5.3 m at the latitude of the scene centre. A degree of longitude shrinks towards the poles while a degree of latitude doesn't, so the pixels are non-square in degrees and the values change from scene to scene. At 40°N, for example, 5.3 m translates to roughly 0.0000621° longitude by 0.0000477° latitude. ### L2A geometry L2A is reprojected into the UTM zone that contains the scene, so pixels are square and measured in metres. Read the zone from the STAC item's `proj:epsg` property rather than assuming one, because scenes in the catalog span many zones. Each item also reports its ground sample distance in `gsd`; values sit close to 5.3 m and vary slightly with off-nadir angle. :::info Full specifications for both levels are in the [Product Guide](https://guide.wyvern.space/data-product-guide/technical-specifications/processing-levels/). ::: ## Working with L1B L1B is the right starting point if you want to run your own atmospheric correction or study the atmosphere itself. For most surface analysis, start from L2A. If you do use L1B, remember what radiance contains. On top of practical steps such as reprojecting to your project's coordinate reference system or clipping to your area of interest, the values still include the shape of the solar spectrum and the effect of the atmosphere. Ozone, water vapour, oxygen, and carbon dioxide absorb light at specific wavelengths, so less of it reaches the sensor there and the signal-to-noise ratio (SNR) in those bands is lower. The ranges where the atmosphere lets light through are called [atmospheric windows](https://gisgeography.com/atmospheric-window/). A vegetation spectrum as L1B radiance (left) and as L2A reflectance after correction with ENVI QUAC (right). The dip near 760 nm on the left is oxygen absorption in the atmosphere. After correction it's gone, and the red edge between 680 and 750 nm is clear. To get from L1B to surface reflectance you need an atmospheric correction tool. The [atmospheric correction](atmospheric_correction) page lists the ones that work with Wyvern data. ## Further reading * [Data Processing Levels (NASA)](https://www.earthdata.nasa.gov/learn/earth-observation-data-basics/data-processing-levels) --- [Source: https://knowledge.wyvern.space/docs/tutorials/python/water_quality_indices] # Water Quality Monitoring with Python ## Exploring Lake Maracaibo, Venezuela Lake Maracaibo, Venezuela is an ecologically stressed water body in South America. It's impacted by algae blooms, turbidity, and runoff from decades of oil extraction and urban discharge. This makes it a great place to see what water quality monitoring looks like from space. Hyperspectral imagery is well suited to this kind of work. Water quality signals like suspended sediment, chlorophyll, and cyanobacteria each leave their own fingerprint across the spectrum, and a hyperspectral sensor captures enough bands to tell them apart. In this tutorial we'll use a Dragonette-003 hyperspectral image captured over Lake Maracaibo on January 4, 2026 to monitor what is happening at the water's surface. By the end you'll know how to: - Load a Wyvern hyperspectral GeoTIFF - Automatically select spectral bands by wavelength - Calculate water quality indices like NDWI, NDTI, and NDCI to map conditions across the scene - Visualize and interpret the results to tell sediment, algae, and clear water apart - Export your findings as a GeoJSON file that can be loaded into QGIS or shared ## Calculating Indices Calculating indices like NDWI or NDCI is an extremely common task in geospatial analysis. At their core, index calculations are just math applied to arrays of numbers — and raster images are exactly that. Each band in a hyperspectral image is a 2D array of reflectance values, one per pixel. Calculating an index means combining two or more of those arrays using a simple formula, chosen so the result highlights something specific in the water like sediment or algae. Many tools like QGIS provide a user-friendly interface for this, but with Python we have a lot more flexibility and control. Here's what we'll do: - Load the image and prepare it for analysis - Select the right spectral bands for each index - Calculate three water quality indices - Visualize and interpret the results ## Packages We'll use the following packages in this tutorial: - **NumPy** is fantastic at array calculations, but it can't read our raster files into arrays on its own. For that, we'll use **Rasterio** (which relies on GDAL), a Python package for reading and writing geospatial raster data that lets us open GeoTIFF files, read band data into NumPy arrays, and access metadata like wavelengths and nodata values. - **Matplotlib** for visualizing results. - **JSON** (built-in) for reading our metadata file. - **GeoPandas** extends **pandas** (Python's standard library for tabular data) to support geographic data. We'll use it first to load our water mask polygon drawn in QGIS, and again, at the end of the tutorial to save our results as a GeoJSON file that can be loaded into QGIS or any other GIS tool. You can install any missing packages with `pip install rasterio numpy matplotlib geopandas`. Try it yourself Runs free in your browser with no install or setup, and pulls the imagery straight from our Open Data Program. Or view the notebook source on GitHub Open in Colab ## Getting started The notebook is self-contained and installs its own dependencies, so there is nothing to set up in advance. - **Google Colab (easiest):** open the notebook using the "Open in Colab" badge above, then run the cells top to bottom. - **Local Jupyter:** open the notebook in JupyterLab or VS Code (Python 3.10+) and run it. Download the notebook and water mask to follow along [here](https://github.com/Nrevyw/wyvern-public-resources/tree/main/tutorial-notebooks/water-quality-indices-maracaibo)! The GeoTIFF, STAC JSON, and data mask needed for this tutorial will be downloaded automatically into the same folder once you run the second cell! ```python import json import numpy as np import rasterio import matplotlib.pyplot as plt import geopandas as gpd ``` ## Download the Tutorial Files Let's download a few more files we'll be working with for this tutorial! We'll be using a hyperspectral GeoTIFF image of our Lake Maracaibo scene, a STAC JSON metadata file, and a data mask. The GeoTIFF contains the image data, the STAC JSON gives us useful metadata about the scene, and the data mask helps identify which pixels contain valid data. Run the cell below to download these files. Once the download is complete, they'll appear in the same folder as your notebook! ```python import requests from pathlib import Path IMAGE_URL = "https://wyvern-data.com/wyvern_dragonette-003_20260104T142819_394db198_l2a/wyvern_dragonette-003_20260104T142819_394db198_l2a.tiff" JSON_URL = "https://wyvern-odp.com/wyvern_dragonette-003_20260104T142819_394db198_l2a/wyvern_dragonette-003_20260104T142819_394db198_l2a.json" DATA_MASK_URL = "https://wyvern-data.com/wyvern_dragonette-003_20260104T142819_394db198_l2a/wyvern_dragonette-003_20260104T142819_394db198_l2a_data_mask.tiff" IMAGE_PATH = Path("wyvern_dragonette-003_20260104T142819_394db198_l2a.tiff") JSON_PATH = Path("wyvern_dragonette-003_20260104T142819_394db198_l2a.json") DATA_MASK_PATH = Path("wyvern_dragonette-003_20260104T142819_394db198_l2a_data_mask.tiff") downloads = [ (IMAGE_URL, IMAGE_PATH), (JSON_URL, JSON_PATH), (DATA_MASK_URL, DATA_MASK_PATH), ] for url, path in downloads: if path.exists(): print(f"{path.name} File already here! Skipping download.") else: print(f"Downloading {path.name}...") response = requests.get(url, headers={"User-Agent": "Mozilla/5.0"}) response.raise_for_status() path.write_bytes(response.content) print("Downloads complete!") ``` ### Loading the Data We'll start by reading the scale factor from the metadata JSON file we downloaded earlier. Wyvern's L2A products represent surface reflectance, but the values are stored as scaled 16-bit integers to reduce file size. Applying the scale factor of `0.0001` converts the stored integer representation into floating-point reflectance values, ranging from 0 to 1. This will make the thresholds, plots, and index calculations we use for water quality easier to interpret. ```python # Open the JSON metadata file and extract the scale factor with open(JSON_PATH) as json_file: metadata = json.load(json_file) # The scale factor converts stored uint16 reflectance values into their floating-point representation scale_factor = metadata["assets"]["Cloud optimized GeoTiff"]["raster:bands"][0]["scale"] # if you don't have the JSON file, you can hardcode the scale factor: # scale_factor = 0.0001 print(f"Scale factor: {scale_factor}") ``` Now let's open the GeoTIFF. We'll replace any nodata pixels with `np.nan` so they are excluded from our calculations, then apply the scale factor to recover the floating-point reflectance values. ```python with rasterio.open(IMAGE_PATH) as image_file: # Convert to float so we can use np.nan and decimal values image_array = image_file.read().astype(np.float32) nodata = image_file.nodata # Replace nodata pixels with NaN so they are excluded from calculations if nodata is not None: image_array[image_array == nodata] = np.nan # Apply the scale factor to recover floating-point reflectance values image_array *= scale_factor print(image_array.shape) ``` ## Building the Wavelength Map Unlike a standard RGB image, each band in a hyperspectral GeoTIFF corresponds to a specific wavelength. We can find this information in the `eo:bands` field under the GeoTIFF asset in the STAC JSON we downloaded earlier, which lists each band's center wavelength. We'll read these out and build a map from wavelength to band index so we can select bands by wavelength rather than hardcoding band numbers. ```python with open(JSON_PATH) as json_file: metadata = json.load(json_file) wavelength_map = {} for i, band in enumerate(metadata["assets"]["Cloud optimized GeoTiff"]["eo:bands"]): wavelength = round(band["center_wavelength"] * 1000) # convert from μm to nm wavelength_map[wavelength] = i print(wavelength_map) ``` ``` {445: 0, 464: 1, 480: 2, 490: 3, 503: 4, 510: 5, 519: 6, 534: 7, 550: 8, 569: 9, 585: 10, 600: 11, 614: 12, 634: 13, 650: 14, 659: 15, 669: 16, 679: 17, 689: 18, 700: 19, 712: 20, 722: 21, 734: 22, 749: 23, 764: 24, 781: 25, 799: 26, 814: 27, 832: 28, 849: 29, 869: 30} ``` ### Selecting Bands by Wavelength Now we'll write a small helper function that finds the nearest available band to any wavelength we ask for. This means if we request 550nm and the closest available band is 549nm, it'll find it automatically. We'll use this in a moment to pull the exact bands each water quality index needs. ```python def get_band(wavelength_nm: int) -> np.ndarray: """ Returns the image band closest to the requested wavelength. Args: wavelength_nm: Target wavelength in nanometres. Returns: A 2D array containing reflectance values for the closest available wavelength band. """ available_wavelengths = list(wavelength_map.keys()) closest_wavelength = min(available_wavelengths, key=lambda w: abs(w - wavelength_nm)) print(f"Requested {wavelength_nm}nm -- using {closest_wavelength}nm") return image_array[wavelength_map[closest_wavelength]] ``` ## Exploring the Scene Before we calculate any indices, let's take a look at our water scene. We'll render it two ways directly from the GeoTIFF using our `get_band()` function to select the right wavelengths: - RGB for a natural colour view - CIR (Colour Infrared) which shifts the bands to highlight vegetation and surface water patterns. For RGB we use the red (~650nm), green (~550nm), and blue (~450nm) bands. These are the same three colour channels your eyes see. For CIR we remove the blue band and add NIR (~800nm), shifting everything along so NIR takes the red channel, red takes green, and green takes blue. This makes vegetation appear bright red since plants reflect NIR strongly. We also normalize each composite before displaying it. Raw reflectance values can have a few very bright or very dark outlier pixels that would make the rest of the image look washed out or too dark. To fix this, we stretch the values between the 2nd and 98th percentile, clipping the brightest 2% and darkest 2% of pixels so the contrast across the rest of the image looks natural. This is only for visualization and doesn't affect any of our calculations. Here, we will normalize each band individually rather than the whole composite at once. This stops a single dominant band (like NIR, which water absorbs very strongly) from pulling down the brightness of the other channels and making the image too dark. ```python # Get bands for RGB and CIR red_band = get_band(650) green_band = get_band(550) blue_band = get_band(450) nir_band = get_band(800) def normalize(array: np.ndarray) -> np.ndarray: """ Stretch an image band to the 0–1 range for display. Values below the 2nd percentile and above the 98th percentile are clipped to reduce the influence of unusually dark or bright pixels. Args: array : np.ndarray A 2D image band containing reflectance values. Returns: np.ndarray The normalized image band with values between 0 and 1. """ vmin = np.nanpercentile(array, 2) vmax = np.nanpercentile(array, 98) return np.clip((array - vmin) / (vmax - vmin), 0, 1) # Normalize each band individually and stack into composites rgb_display = np.dstack([normalize(red_band), normalize(green_band), normalize(blue_band)]) cir_display = np.dstack([normalize(nir_band), normalize(red_band), normalize(green_band)]) # Plot side by side fig, axes = plt.subplots(1, 2, figsize=(14, 8)) axes[0].imshow(rgb_display) axes[0].set_title('RGB', fontsize=12, fontweight='bold') axes[0].axis('off') axes[1].imshow(cir_display) axes[1].set_title('CIR (Colour Infrared)', fontsize=12, fontweight='bold') axes[1].axis('off') plt.tight_layout() plt.show() ``` ![RGB and Colour Infrared composites of the Lake Maracaibo scene, side by side](images/maracaibo_scene_overview.png 'RGB and CIR composites of Lake Maracaibo :size=900') The RGB image gives us a natural colour view of the scene. The strait runs through the centre, with the city of Maracaibo visible on the left. The bright green swirling patterns in the water suggest some kind of biological or suspended material activity, possibly algae-related. In the CIR image, vegetation on land appears bright red as expected. The algae-like features in the water show up as bright pink/magenta streaks. This is because whatever is causing those patterns appears to be reflecting NIR light, which is unusual for typical water and may suggest the presence of biological material. These images give us a good starting point but can only tell us so much. To really dig into what's happening in the water, we'll use hyperspectral indices to map specific signals like turbidity and chlorophyll across the scene. Let's get started! ## The Indices Here's an overview of each index we'll calculate and what it tells us about the water. When choosing indices for a scene, a good starting point is to think about what you're trying to measure and what's physically happening in the water. For Lake Maracaibo, where algae blooms and turbidity are both known to be present, these indices can give us a clearer picture of the water's condition. You can find a list of indices in the [Wyvern Hyperspectral Index Library](/hyperspectral_library). **NDWI (Normalized Difference Water Index)** `(green - nir) / (green + nir)` (green: 550nm, NIR: 800nm) Uses the green and NIR bands to identify water. Water commonly has higher NDWI values than surrounding land because it absorbs strongly in the near-infrared. However, the exact values vary with scene conditions, atmospheric effects, turbidity, algae, and the selected bands. It's often used as a water mask before applying other water quality indices. **NDTI (Normalized Difference Turbidity Index)** `(red - green) / (red + green)` (red: 650nm, green: 550nm) Uses the red and green bands to estimate turbidity. Turbid water carrying suspended sediment scatters more red light relative to green, so higher NDTI values indicate more turbid water. **NDCI (Normalized Difference Chlorophyll Index)** `(red_edge - red) / (red_edge + red)` (red-edge: 712nm, red: 650nm) Uses the red-edge and red bands to detect chlorophyll-a in water. This is where hyperspectral imagery really shines. Dragonette's narrow bands give us coverage right at the wavelength where the chlorophyll signal is strongest. ### Getting our Bands Now we can use our helper function to select the right bands for each water quality index. We just pass in the wavelength we need and it finds the closest available band automatically. ```python green = get_band(550) # NDWI, NDTI nir = get_band(800) # NDWI red = get_band(650) # NDTI, NDCI re = get_band(712) # NDCI ``` ### Calculating the Indices Now we can calculate each index using NumPy. Each equation is applied across every pixel in the image at once — so no loops needed! ```python # Calculate each index across every pixel ndwi = (green - nir) / (green + nir) ndti = (red - green) / (red + green) ndci = (re - red) / (re + red) ``` ## NDWI: Mapping Water Let's start by plotting NDWI across the full scene. This gives us a clear picture of where the water is — the first step before we start assessing what's in it. Rather than hardcoding `vmin` and `vmax`, we'll use the 2nd and 98th percentiles of the data. This automatically clips out extreme outlier values at both ends and makes the variation across the scene much more visible. ```python ndwi_vmin = np.nanpercentile(ndwi, 2) ndwi_vmax = np.nanpercentile(ndwi, 98) print(f"NDWI minimum: {np.nanmin(ndwi):.4f}") print(f"NDWI maximum: {np.nanmax(ndwi):.4f}") print(f"NDWI 2nd percentile: {ndwi_vmin:.4f}") print(f"NDWI 98th percentile: {ndwi_vmax:.4f}") plt.figure(figsize=(5, 8.5)) plt.imshow(ndwi, cmap='RdYlBu', vmin=ndwi_vmin, vmax=ndwi_vmax) plt.colorbar(label='NDWI') plt.title('NDWI') plt.axis('off') plt.show() ``` ![NDWI map of the Lake Maracaibo scene](images/maracaibo_ndwi.png 'NDWI across the full scene :size=500') Blue areas are water, red areas are land. Within the water body you can see some areas showing lower NDWI values than the open water, possibly where surface conditions like algae are affecting the green and NIR reflectance differently from clear water. ### Water Mask NDWI is commonly used as a water mask, but it struggles with this scene. Algae reflects both green and NIR light, which can push NDWI values lower than expected for water, or even become negative in heavily algae-covered areas. Using a simple NDWI threshold risks masking out exactly the water pixels we're most interested in and the algae-affected water is the whole reason we're here. Instead, we'll use two masks combined: - **A water polygon** manually drawn in QGIS following the shoreline of the scene, which cleanly separates water from land regardless of what's happening spectrally in the water. - **Wyvern's QA_CLEAR_MASK**, a band included with every L2A product that flags whether each pixel is clear of cloud, haze, and shadow. A value of `1` means the pixel is clear and good to use, `0` means it's affected by atmospheric interference. This automatically excludes the cloudy pixels visible at the bottom of our scene. Together these will give us a clean mask of clear, cloud-free water pixels. ```python from rasterio.features import geometry_mask WATER_MASK_PATH = "water_mask.geojson" # Load the water polygon drawn in QGIS, and extract geometries for raster masking water_mask = gpd.read_file(WATER_MASK_PATH) geometries = water_mask.geometry # Load Wyvern's QA_CLEAR_MASK: 1 = clear pixel, 0 = cloud/haze/shadow with rasterio.open(DATA_MASK_PATH) as data_mask_file: qa_clear = data_mask_file.read(1) # Rasterize the water polygon to match the image dimensions with rasterio.open(IMAGE_PATH) as image_file: water_polygon_mask = geometry_mask( geometries, transform=image_file.transform, invert=True, out_shape=(image_file.height, image_file.width) ) # Combine water polygon and QA clear mask and keep only clear water pixels water_mask = water_polygon_mask & (qa_clear == 1) ``` Let's visualize the water mask to make sure it looks correct before applying it to our water quality indices. ```python plt.figure(figsize=(5, 8.5)) plt.imshow(water_mask, cmap='Blues') plt.title('Water Mask') plt.axis('off') plt.show() ``` ![Binary water mask of the Lake Maracaibo scene, combining a manually drawn polygon with Wyvern's QA_CLEAR_MASK](images/maracaibo_water_mask.png 'Water mask :size=500') ## NDTI: Turbidity NDTI uses the red and green bands to estimate turbidity. Higher values indicate more suspended sediment in the water column, while lower values suggest less turbidity. We calculated it earlier using: `ndti = (red - green) / (red + green)` For this scene we'd expect to see higher turbidity near the shoreline and in shallower areas where sediment is stirred up, and lower values in the deeper open water of the strait. Let's check the value range before plotting. This helps us set `vmin` and `vmax` so the variation in the water is visible and doesn't get lost in the full -1 to 1 range. ```python # Apply the water mask to NDTI, excludes land pixels from the plot ndti_masked = np.where(water_mask, ndti, np.nan) # Check the value range to set colormap limits print(f"Min: {np.nanmin(ndti_masked):.4f}") print(f"Max: {np.nanmax(ndti_masked):.4f}") print(f"2nd percentile: {np.nanpercentile(ndti_masked, 2):.4f}") print(f"98th percentile: {np.nanpercentile(ndti_masked, 98):.4f}") ``` ```python # Use percentiles to take out extreme values and improve contrast vmin = np.nanpercentile(ndti_masked, 2) vmax = np.nanpercentile(ndti_masked, 98) plt.figure(figsize=(5, 8.5)) plt.imshow(ndti_masked, cmap='BrBG_r', vmin=vmin, vmax=vmax) plt.colorbar(label='NDTI') plt.title('NDTI') plt.axis('off') plt.show() ``` ![NDTI map of the Lake Maracaibo scene, masked to water pixels only](images/maracaibo_ndti.png 'NDTI across the water mask :size=500') The NDTI map shows some really clear spatial patterns! Higher values (brown and tan) near the shorelines may reflect shallower areas close to the coast, not necessarily more turbid water. However, the lower values (teal and green-blue) running through the centre of the strait are interesting. NDTI can't tell the difference between turbidity and other surface conditions, so those lower values might not just mean cleaner, less turbid water — there could be something else going on spectrally. Based on what we've seen from the RGB visual, those lower values could reflect the influence of algae or other biological material. ### Highlighting Areas of High Turbidity The original NDTI map shows the full range of turbidity values across the scene. However, this can make the most turbid areas harder to identify because lower and moderate values still occupy much of the colour scale. To make the strongest turbidity patterns stand out, we'll display only the highest 20% of NDTI values within the water area. ```python # Only show the highest NDTI values so high-turbidity areas stand out threshold = np.nanpercentile(ndti_masked, 80) # Hide pixels below the threshold high_ndti = np.where(ndti_masked >= threshold, ndti_masked, np.nan) # Scale the colour range using only the highlighted values vmin = threshold vmax = np.nanpercentile(ndti_masked, 98) plt.figure(figsize=(5, 8.5)) plt.imshow(high_ndti, cmap='YlOrRd', vmin=vmin, vmax=vmax) plt.colorbar(label='NDTI') plt.title('High Turbidity Areas') plt.axis('off') plt.show() ``` ![NDTI map showing only the highest 20% of turbidity values within the water mask](images/maracaibo_high_ndti.png 'Highlighted high-turbidity areas :size=500') Scaling the colormap this way lets the high-value areas stand out rather than blending into a gradient across the scene. The strongest signal shows up along the shorelines and in the shallower areas close to the coast, though as we noted, this doesn't necessarily mean more turbid water. Let's focus on the lower signal in the open channel running through the centre. Something other than sediment might be influencing those pixels — algae absorbs red light and reflects green, which can push NDTI low in areas with significant biological surface cover. NDTI alone can't tell the difference. We'll look into this further with NDCI! ## NDCI: Chlorophyll NDCI uses the red-edge and red bands to detect chlorophyll-a in water. We calculated it earlier using: `ndci = (re - red) / (re + red)` Hyperspectral imagery can be really helpful here because it can give us narrow bands across the visible and red-edge regions! Chlorophyll absorbs strongly in the blue and red parts of the spectrum, around ~440 nm and ~670 nm. NDCI uses the contrast between the red absorption region and the red-edge region around ~700–720 nm, where reflectance increases due to scattering from algal and cellular structures. Higher NDCI values means a stronger chlorophyll signal in the water. For this scene, if the lower NDTI values we saw in the channel are actually caused by algae, we'd expect to see higher NDCI values in those same areas. Let's check the value range before plotting: ```python # Apply the water mask to NDCI ndci_masked = np.where(water_mask, ndci, np.nan) # Check the value range to set colormap limits print(f"Min: {np.nanmin(ndci_masked):.4f}") print(f"Max: {np.nanmax(ndci_masked):.4f}") print(f"2nd percentile: {np.nanpercentile(ndci_masked, 2):.4f}") print(f"98th percentile: {np.nanpercentile(ndci_masked, 98):.4f}") ``` ```python # Use percentiles to clip extreme values and improve contrast vmin = np.nanpercentile(ndci_masked, 2) vmax = np.nanpercentile(ndci_masked, 98) plt.figure(figsize=(5, 8.5)) plt.imshow(ndci_masked, cmap='RdYlGn', vmin=vmin, vmax=vmax) plt.colorbar(label='NDCI') plt.title('NDCI') plt.axis('off') plt.show() ``` ![NDCI map of the Lake Maracaibo scene, masked to water pixels only](images/maracaibo_ndci.png 'NDCI across the water mask :size=500') The NDCI map gives us a clearer picture. The bright green areas running through the strait line up exactly with the swirling patterns we saw in the RGB image. Higher NDCI values here mean elevated chlorophyll, and could point to increased algal activity. This also helps make sense of what we saw in the NDTI map. Those lower turbidity values in the channel might not just mean cleaner water. Algae and other biological material can affect the red/green ratio too. The open water in the upper portion of the scene shows lower NDCI values (red and orange), suggesting less chlorophyll which makes sense for more open, less nutrient-rich water. This is a good example of why using multiple indices together gives you a much more complete picture of the water's condition than any single index on its own! ## Areas of Interest Let's zoom into two specific areas to see what's going on in the water more closely: one where turbidity is the dominant signal, and one where chlorophyll is. For each area we'll look at four views side by side: RGB and CIR for visual context, then NDTI and NDCI to see how the indices respond. Comparing them directly is a good way to build intuition for what each index is actually picking up and where they tell different stories about the water. ![RGB, CIR, NDTI, and NDCI comparison panels for the high-turbidity area](images/maracaibo_turbidity_comparison.png 'High-turbidity area comparison :size=900') ![RGB, CIR, NDTI, and NDCI comparison panels for the high-chlorophyll area](images/maracaibo_chlorophyll_comparison.png 'High-chlorophyll area comparison :size=900') ### High-Turbidity Area In the RGB image, there's a faint plume visible in the water that's quite easy to miss at first. The CIR adds a little more context, with a slightly different tone in that region hinting that something is going on spectrally, though it's still not obvious what. NDTI is where it becomes clear. The plume lights up with a strong turbidity signal, confirming that whatever is in the water is scattering red light, which is consistent with suspended sediment. And when we check NDCI, chlorophyll is low across the whole area. So this isn't algae. It could be that something has physically disturbed the water here and stirred up sediment from below. A passing ship is a reasonable explanation, and there is evidence of possible vessel traffic in the RGB if you look closely. ### High-Chlorophyll Area Here the RGB tells a much more obvious story straight away. The bright green swirling streaks are hard to miss. In CIR they show up as magenta, which tells us the surface is reflecting strongly in NIR. That's unusual for water and points toward biological material at the surface. NDTI is strongly negative across the whole area, the opposite of the turbidity region. Low red reflectance relative to green is consistent with chlorophyll absorbing red light. And NDCI confirms it: the swirling patterns map almost perfectly onto high chlorophyll values. This is very good evidence of an algal bloom. Together, these two areas are a great example of how combining indices lets us tell one water quality concern from another — sediment in one spot, an algal bloom in the other. ## Comparing Spectral Curves So far, we've been looking at the scene through image displays and index maps. Another useful way to explore water quality with hyperspectral data is to look at the spectra directly. Each condition in the water has its own shape across wavelength, and plotting the curve makes that fingerprint easy to see. In this section, we'll compare the reflectance curves from three points in the water: one in a high-turbidity area, one in a high-chlorophyll area, and one in clear open water. We'll sample a single representative pixel for each condition and plot its curve across Wyvern's 31 bands. The image below shows where each sample point sits within the scene. This helps us connect the maps back to the underlying data. Instead of only seeing where an index is high or low, we can look at how each condition behaves across wavelength and see which parts of the spectrum are driving the patterns we observed. ![Sample pixels for the high-turbidity, high-chlorophyll, and open water spectral curves, shown on the RGB scene](images/maracaibo_sample_pixels.png 'Sample pixels for spectral curves :size=600') ![Spectral reflectance curves for the high-turbidity, high-chlorophyll, and open water sample pixels across all 31 bands](images/maracaibo_spectral_curves.png 'Spectral curves from selected pixels :size=900') The high-chlorophyll pixel (green) peaks in the green band around 550nm, dips through the red where chlorophyll absorbs, then rises sharply starting around 700 nm, climbing far above the other two and staying high all the way through the near-infrared. That steep jump is the chlorophyll red-edge: chlorophyll absorbs red light but reflects strongly just past it, and it's the feature NDCI is built to detect. A strong red-edge like this points to a lot of chlorophyll material at the surface and is the signal you'd expect from dense vegetation or an algae-rich patch of water. The high-turbidity pixel (brown) sits above open water through most of the visible range. That makes sense since suspended sediment scatters light broadly across the spectrum rather than at one specific wavelength, so turbid water is brighter overall in the visible. It doesn't show the red-edge jump: in the near-infrared it stays low and flat, close to open water. That contrast is the difference between a turbidity signal and a chlorophyll signal, and it's why NDTI looks at the red vs green balance while NDCI looks at the red-edge. Open water (blue) is the lowest curve throughout and gives us a baseline. It stays low across the visible and has none of the red-edge rise we see in the bloom. Clear water absorbs more of the light that hits it, while the suspended sediment in turbid water scatters light more. These spectra confirm that what NDTI and NDCI flagged really are different things happening in the water! ## Sharing Your Results So far we've visualized the indices as maps inside the notebook. But if you want to share the results with someone who doesn't have Python, load them into QGIS, or feed them into another workflow, we can export them as a GeoJSON file. Here we'll do that with our NDCI results, turning the areas with the highest chlorophyll signal into a vector layer of polygons that can be used anywhere. We'll treat the top 10% of NDCI values within the water mask as high-chlorophyll pixels, then use `rasterio.features.shapes()` to trace their edges into geographic polygons and save them as a GeoJSON using GeoPandas. ```python from rasterio.features import shapes from shapely.geometry import shape # Select the top 10% of NDCI values as high chlorophyll threshold = np.nanpercentile(ndci_masked, 90) high_chlorophyll_mask = np.where(ndci_masked >= threshold, 1, 0).astype(np.uint8) # Get the image transform and CRS with rasterio.open(IMAGE_PATH) as image_file: transform = image_file.transform crs = image_file.crs # Trace pixels into polygons high_chlorophyll_polygons = [ shape(polygon_geometry) for polygon_geometry, pixel_value in shapes(high_chlorophyll_mask, transform=transform) if pixel_value == 1 ] # Save as GeoJSON in WGS84 so it loads correctly in QGIS and other GIS tools gdf = gpd.GeoDataFrame(geometry=high_chlorophyll_polygons, crs=crs) gdf.to_crs("EPSG:4326").to_file("maracaibo_high_chlorophyll.geojson", driver="GeoJSON") print("Saved maracaibo_high_chlorophyll.geojson!") ``` ``` Saved maracaibo_high_chlorophyll.geojson! ``` The GeoJSON file can now be dragged straight into QGIS and overlaid on a satellite basemap. The polygons trace the high-chlorophyll areas we identified from the NDCI map. ![High chlorophyll areas exported as GeoJSON and loaded into QGIS over a Google Satellite basemap.](images/maracaibo_high_chlorophyll_qgis.png 'High chlorophyll areas in QGIS :size=900') *High chlorophyll areas exported as GeoJSON and loaded into QGIS over a Google Satellite basemap.* ## Other Packages NumPy is the most straightforward tool for index calculations and what we've used throughout this tutorial, but there are several other Python packages worth knowing about! ### HyperCoast [HyperCoast](https://hypercoast.org/) is a Python package for visualizing and analyzing hyperspectral data, with a focus on coastal and water applications — so it's a natural fit for water-quality work like this. It lets you view hyperspectral scenes, swap band combinations, pull out spectral signatures, and even explore your data in 3D. There's also a QGIS plugin if you'd rather work interactively with less code. Here's a [walkthrough for visualizing Wyvern open data](https://www.youtube.com/watch?v=BYrQdhz5k5o) to get you started. It's peer-reviewed through JOSS, and you can install it with `pip install hypercoast`. ### NumExpr [NumExpr](https://numexpr.readthedocs.io/en/latest/) speeds up NumPy array operations by evaluating expressions in chunks and taking advantage of multiple CPU cores. For large images it can be significantly faster than plain NumPy. An NDWI calculation that would look like `(green - nir) / (green + nir)` in NumPy becomes `ne.evaluate("(green - nir) / (green + nir)")` in NumExpr — the equation is the same, just passed as a string. ### EarthPy [EarthPy](https://earthpy.readthedocs.io/en/latest/api/earthpy.spatial.html#earthpy.spatial.normalized_diff) provides a `normalized_diff()` function which can be used to calculate simple two-band indices like NDWI and NDCI without writing the equation yourself. ### Xarray-Spatial [Xarray-Spatial](https://xarray-spatial.org/user_guide/multispectral.html) has built-in functions for calculating a number of useful indices using Xarray arrays, which can be handy for working with larger datasets. ### Awesome Spectral Indices (Spyndex) David Montero Loaiza's [Awesome Spectral Indices](https://github.com/awesome-spectral-indices/awesome-spectral-indices) project contains a standardised database of spectral indices. For Python, [Spyndex](https://github.com/awesome-spectral-indices/spyndex) lets you pull equations directly from the database rather than writing them yourself — worth exploring if you plan to calculate many different indices regularly. ## What's Next? In this tutorial we used three indices to explore water quality patterns in a single scene. But this is just a starting point — there are many more indices in the [Wyvern Hyperspectral Index Library](/hyperspectral_library) worth exploring. A few ideas to take this further: - Try applying these indices to a different Wyvern scene and compare the results. Browse available scenes on the [Wyvern Open Data Program](https://opendata.wyvern.space) - Look into the [Visualizing Spectra & Plotting](/docs/tutorials/python/visualizing_spectra) tutorials for more ways to explore hyperspectral data in Python --- [Source: https://knowledge.wyvern.space/docs/tutorials/qgis/creating_indices] # Index Calculations in QGIS Natural colour RENDVI Wilderness area south of Halifax, Canada: Natural colour vs. RENDVI ### Before begining This tutorial assumes the following steps have already been taken: - QGIS has been downloaded and installed on your machine. Visit [qgis.org](https://www.qgis.org/) to download QGIS. - You have loaded a Wyvern image into QGIS, feel free the explore the [open data program](https://opendata.wyvern.space/#/?.language=en) for one that interests you. - We will be using [this image](https://opendata.wyvern.space/#/surface-reflectance/wyvern_dragonette-003_20250912T143415_6ead8ea3_l2a/wyvern_dragonette-003_20250912T143415_6ead8ea3_l2a.json?.language=en) over southern Halifax, Canada for this tutorial. --- QGIS provides a nice user interface for many workflows GIS analysts will need. One of the most common workflows when working with satellite imagery is creating indices. Calculating an index in QGIS is as simple as using the built in raster calculator to compute pixel values using selected bands based on the index you want to produce! ### Indices Data scientists and analysts create indices from satellite imagery to simplify complex spectral data into clear, interpretable measures of surface conditions. Indices such as the Red Edge Normalized Difference Vegetation Index (RENDVI) combine multiple spectral bands to reveal patterns in vegetation health that are not obvious in raw imagery. This simplification improves comparability across space and time and supports efficient large-scale monitoring and analysis. For this tutorial we will compute the RENDVI over a wilderness area south of Halifax on the east coast of Canada. RENDVI is just one of many index calculations enabled by hyperspectral data, we have tons of examples of indices available for you to explore at our [hyperspectral library](/hyperspectral_library). ### The Workflow To start, we will open our Dragonette data in QGIS and lightly style it to give us a natural color composite. Natural color composite of wilderness area south of Halifax (Bands: 11, 5, 1) :::info Note: Your image may look different on load. If you are having trouble loading in or creating an RGB image, please follow the steps provided in [this tutorial](loading_data).\ Depending on your machine, it may make sense to [clip the image](https://docs.qgis.org/3.40/en/docs/user_manual/processing_algs/gdal/rasterextraction.html#clip-raster-by-extent) to save on memory. ::: Now that we have our image loaded in, we can calculate our RENDVI.\ RENDVI is a fairly simple index calculated using two red bands, one at 750nm ($R_{750}$) and the other at 710nm ($R_{710}$). $$\text{RENDVI} = \frac{R_{750} - R_{710}}{R_{750} + R_{710}}$$ Red Edge Normalized Difference Vegetation Index (RENDVI) Equation We can refer to the [RENDVI entry in our hyperspectral library](/hyperspectral_library#red-edge-ndvi-rendvi) and see that we should use `Band 24` for our $R_{750}$ band and `Band 21` for our $R_{710}$ band. Translating this to QGIS means opening the `Raster Calculator` (`Raster -> Raster Calculator`). We then see `Dragonette-003`'s 31 bands that we can choose from. All we have to do now is plug in our `Band 24` as the $R_{750}$ and our `Band 21` as the $R_{710}$!\ The formula in your raster calculator should look something like this: ``` ("wyvern_dragonette-003_20250912T143415_6ead8ea3_l2a@24" - "wyvern_dragonette-003_20250912T143415_6ead8ea3_l2a@21") / ("wyvern_dragonette-003_20250912T143415_6ead8ea3_l2a@24" + "wyvern_dragonette-003_20250912T143415_6ead8ea3_l2a@21") ``` Next just name our `Output layer` something sensible, choose your prefered `Output format` (GeoTIFF works well), and click ok! `Raster Calculator` will run that formula on each pixel in our image, and create a new raster with values ranging from -1 to 1. `Raster Calculator` with RENDVI calculation plugged in Now we have our RENDVI raster calculated and displayed in our workspace! We can apply some symbology to our layer to be more consistent with typical RENDVI visualizations. To do so: 1. Double click the newly created raster in the `layers` panel. 2. In the `symbology` tab, change the `Render type` to `Singleband pseudocolor`. 3. Here we will choose our color ramp, RENDVI typically uses a red-yellow-green color ramp. - You can find this by clicking the down arrow on the right of `Color Ramp` - `Color ramp -> All Color Ramps -> RdYlGn` Applying symbology to our RENDVI layer :::tip Pro Tip You can determine your own min/max cutoffs for styling by computing a histogram of the raster values by going to\ `Histogram -> Compute Histogram`.\ We will use `Min: 0, Max: 0.4` for our RENDVI visualization. ::: Voila! You have now created an RENDVI layer out of Wyvern's hyperspectral data! You'll notice areas that pixels containing any water or synthetic materials in the image will appear red, while healthier vegetation will appear darker green. Stylized RENDVI output layer Now that you know how to translate formulas for hyperspectral indices into QGIS's raster calculator, we implore you to explore other images in our [open data program](https://opendata.wyvern.space/#/?.language=en) and apply other indices from our [hyperspectral library](/hyperspectral_library). Hyperspectral enables more precise index calculation, allowing for better data for scientists, policy makers, and you. Unlock the full story hidden in the spectrum and start exploring what’s possible. --- ### Further Readings - [Plant Health Indices](../../documentation/datasets_and_use_cases/plant_health_indices) --- [Source: https://knowledge.wyvern.space/docs/documentation/getting_started/stac_documentation_and_tools] # Understanding the STAC format Every Wyvern scene comes with a [SpatioTemporal Asset Catalog (STAC)](https://stacspec.org/en/about/) record, and the Open Data Program is organized as a STAC catalog. This page explains the three building blocks of STAC, shows how our catalog uses them, and lists the tools that read it. Once you know the structure, you can find scenes by place, date, or application without downloading anything. ## What is STAC? STAC is an open standard for describing geospatial data so it can be searched and combined consistently. Because every provider that uses STAC describes its data the same way, the code you write to find Wyvern scenes will also find Sentinel-2 or Landsat scenes with few changes. A *spatiotemporal asset* is any file that records information about the Earth at a particular place and time: satellite, drone, or aerial imagery of any kind, LiDAR point clouds, or video. A STAC catalog is a set of JSON files that link to one another. Each file follows the core specification for the kind of object it describes, and the specification can be extended with extra fields for particular kinds of data `[1]`. Following the links takes you from the top of the catalog down to the assets themselves: in our case, the hyperspectral images, their masks, and their previews. ## Key components of STAC ### Item An Item describes one scene. It records when and where the scene was captured and links to its files, such as the GeoTIFF and thumbnails. An Item is a [GeoJSON](https://en.wikipedia.org/wiki/GeoJSON) feature with extra metadata fields, so any modern GIS or geospatial library can read its footprint directly. > **Example**: [`wyvern_dragonette-003_20260104T142819_394db198_l2a`](https://wyvern-odp.com/wyvern_dragonette-003_20260104T142819_394db198_l2a/wyvern_dragonette-003_20260104T142819_394db198_l2a.json), a Dragonette-003 scene captured on January 4, 2026, used in the [water quality tutorial](/docs/tutorials/python/water_quality_indices). The Item specification uses the standard GeoJSON fields and adds its own. Figure 1 shows them. Figure 1: The STAC Item specification. Fields in grey (stac_extensions and collection) are optional. `[2]` A Wyvern Item links to seven assets: the hyperspectral COG, a Data Mask, a Pixel Quality Mask, an overview image, a thumbnail, the item's own metadata, and a zip file of everything. Its properties record the processing level, platform, cloud cover, and sun and viewing angles, and the COG asset lists every band's centre wavelength and FWHM. ### Catalog A Catalog organizes Items into a hierarchy you can browse, like a table of contents. A `catalog.json` file links to any combination of other Catalogs, Collections, and Items. Figure 2 shows its required fields. > **Example**: The [Open Data root catalog](https://wyvern-odp.com/catalog.json), which groups the same scenes four ways: by year, by application (agriculture, coastal, forestry, mining, and others), by spectral range, and a separate set of L1B scenes. Figure 2: The STAC Catalog specification. Fields in white (title and stac_extensions) are optional. `[2]` ### Collection A Collection is a Catalog with added metadata. It groups Items that share characteristics, such as a sensor, a time period, or an acquisition campaign, and it adds required fields that describe the whole group: spatial and temporal extent, licence, and providers. A Catalog just organizes; a Collection also tells you what's in it. Figure 3 shows the extra fields. > **Example**: The [Extended VNIR collection](https://wyvern-odp.com/product-type/extended/collection.json), which groups every 31-band scene from Dragonette-002 onward. Figure 3: The STAC Collection specification. In addition to the optional Catalog fields (title, stac_extensions), providers, keywords, assets, and summaries are optional for Collections. `[2]` ## How Wyvern delivers data with STAC Each Wyvern scene is a Cloud-Optimized GeoTIFF with its own embedded metadata, plus a STAC Item JSON file alongside it. The GeoTIFF gives you the pixels and the band wavelengths. The STAC Item gives you everything else, including the scale factor that turns L2A's stored integers into reflectance. The GeoTIFF doesn't store that factor, so take it from the STAC Item. See [processing levels](processing_levels) for why that matters. :::info Browse the catalog in the [Open Data Program](https://opendata.wyvern.space/), or [contact us](mailto:info@wyvern.space) about tasking new imagery. ::: ## Tools for working with STAC These are the tools we'd reach for first: * [PySTAC](https://pystac.readthedocs.io/en/stable/): reads, writes, and edits STAC Catalogs, Collections, and Items in Python. Use it to walk a static catalog like ours. * [pystac-client](https://pystac-client.readthedocs.io/en/stable/): searches STAC APIs by area, date, and properties from Python. Use it for providers such as Microsoft Planetary Computer or Earth Search that serve Sentinel-2 and Landsat through an API. * [STAC Browser](https://github.com/radiantearth/stac-browser): a web interface for browsing and searching STAC catalogs. The Open Data Program site is built on it. * [GDAL](common_tools#gdal): not a STAC tool as such, but it reads the COGs a STAC Item points to, directly over HTTP, which makes it easy to drop into a STAC workflow. For a working example in Python, the [agents page](../agents) has a short script that walks our catalog, picks a scene, and resolves bands by wavelength. ### Further reading * [STAC documentation](https://stacspec.org/en/about/) * [STAC tutorials](https://stacspec.org/en/tutorials/index.html) ## References `[1]` STAC API Community Standard Work Item Justification. [https://portal.ogc.org/files/105786](https://portal.ogc.org/files/105786) `[2]` Intro to STAC [https://stacspec.org/en/tutorials/intro-to-stac/](https://stacspec.org/en/tutorials/intro-to-stac/) --- [Source: https://knowledge.wyvern.space/docs/documentation/getting_started/hyperspectral_signatures] # Hyperspectral signatures Every pixel in a hyperspectral image has a spectrum: its reflectance in each band. That spectrum depends on what's in the pixel, and comparing it against known spectra is how you work out what the pixel contains. This page covers spectral signatures, the complication of mixed pixels, and the two broad ways of working with them. ## Spectral signatures A spectral signature is a material's reflectance plotted against wavelength. Signatures are usually measured in a laboratory or the field on a pure sample, and collections of them are called spectral libraries. See [material reflectance libraries](../more/material_reflectance) for where to find them. Signatures of whitebark pine and English oak from the USGS Spectral Library version 7. Both are vegetation and share the same overall shape; the differences between them are what a hyperspectral sensor can pick up and a multispectral one mostly can't. Before you compare a library spectrum with a pixel, resample it to the sensor's bands. Library spectra are often sampled every 1 nm, while each Dragonette band averages over 15 to 32 nm, so the two only line up once the library spectrum has been averaged through each band's response. The [rare earth elements tutorial](/docs/tutorials/python/detecting_rare_earth_elements) shows how. ### Mixed pixels A pixel often contains more than one material, and its spectrum is then a blend of theirs, roughly in proportion to how much of the pixel each one covers. This is called a mixed pixel, and it means a pixel's spectrum won't match any single library signature exactly. Coarse pixels (black boxes) over two fields. Pixels inside a field are pure; pixels over the road, railway, and field edges are mixed. The coarser the pixel, the more of them are mixed, because each one covers more ground. At 5.3 m, a Dragonette pixel is small enough to fall inside a single field, roof, or stand of trees far more often than a 30 m pixel is, but mixing never goes away entirely. Field edges, roads, and small features will always produce mixed pixels. There are two ways to deal with mixing: model it, or classify each pixel by its closest match. ### Unmixing pixels The first approach models the mixing directly. *Endmember extraction* finds the purest spectra in a scene, the endmembers, which represent its main materials. *Spectral unmixing* then estimates what fraction of each pixel is made up of each endmember, so instead of labelling a pixel "field" or "road" you get a map of how much of each it contains. [Spectral Python (SPy)](https://www.spectralpython.net/algorithms.html#endmember-extraction) implements endmember extraction and documents how to apply it. ### Classifying pixels The second approach is to compare each pixel with known signatures and assign it to the closest match. Because library spectra are measured under laboratory conditions and image spectra aren't, you're looking for the most similar signature rather than an exact one. The spectral angle mapper (SAM) is the classic example. It compares the *shape* of two spectra and ignores their overall brightness, so a shaded pixel still matches a sunlit library spectrum of the same material. The [rare earth elements tutorial](/docs/tutorials/python/detecting_rare_earth_elements) uses a related detector, ACE, to find bastnaesite, a neodymium-bearing mineral, at the Mountain Pass mine in California. [Spectral Python](https://www.spectralpython.net/algorithms.html#supervised-classification) also provides supervised classification algorithms. --- [Source: https://knowledge.wyvern.space/docs/tutorials/python/detecting_rare_earth_elements] # Detecting Rare Earth Elements The goal of this tutorial: an Nd detection map (right) built from a Wyvern scene over the Mountain Pass REE mine (left). Try it yourself Runs free in your browser with no install or setup, and pulls the imagery straight from our Open Data Program. Or view the notebook source on GitHub Open in Colab Rare earth elements (REEs) power everything from electric vehicles to wind turbines, and finding new sources is a strategic priority. It turns out several REEs leave a spectral fingerprint we can see from orbit. **Neodymium (Nd)** in particular produces sharp, characteristic absorption features in the visible / near-infrared (VNIR) at roughly **585, 745, 810 and 870 nm**, and — crucially — those positions barely move regardless of the host mineral. That makes Nd a reliable proxy for REE-bearing material. This tutorial is inspired by a great paper: Asadzadeh, Koellner & Chabrillat (2024), [*Detecting rare earth elements using EnMAP hyperspectral satellite data: a case study from Mountain Pass, California*](https://doi.org/10.1038/s41598-024-71395-2). It reports the first direct detection of Nd from a spaceborne hyperspectral satellite (earlier REE mapping had been done from the lab, the ground, and aircraft, but not from orbit). This tutorial will create an approachable and easy to run workflow in three small steps: 1. **Continuum removal** to isolate the absorption features from overall brightness. 2. **Spectral resampling** of a real USGS bastnaesite reference spectrum onto Wyvern's bands. 3. **ACE** (Adaptive Cosine Estimator), a matched-filter-style detector, to map Nd-bearing pixels. :::warning Use L2A surface reflectance REE absorptions are subtle. This workflow needs **L2A surface reflectance** (atmospherically corrected) imagery. If you run it on L1B top-of-atmosphere radiance, the solar spectrum and the oxygen absorption band near 760 nm will swamp the features. All imagery on the Open Data Program is L2A unless otherwise noted. ::: ## Getting started The companion notebook is self-contained and installs its own dependencies, so there is nothing to set up in advance. - **Google Colab (easiest):** open the notebook using the "Open in Colab" badge at the top of it, then run the cells top to bottom. - **Local Jupyter:** open the notebook in JupyterLab or VS Code (Python 3.10+) and run it. Either way, the first cell installs the packages that are not already present: ```python %pip install -q rasterio spectral scipy pyarrow ``` ## The image we're using We use a [Wyvern Dragonette-002 scene over Mountain Pass, California](https://opendata.wyvern.space/#/wyvern_dragonette-002_20250422T221452_2bc40137_l2a/wyvern_dragonette-002_20250422T221452_2bc40137_l2a.json), available as an L2A surface-reflectance product on the [Wyvern Open Data Program](https://opendata.wyvern.space/#/?.language=en). The notebook downloads it automatically from its STAC item. ## The Nd fingerprint Neodymium in its common +3 oxidation state (Nd(III)) has a set of narrow absorption bands whose centre wavelengths barely shift with the host mineral, which is what makes them such a dependable fingerprint. Four fall within Wyvern's Extended VNIR range: | Feature | Wavelength | Notes | | ------- | ---------- | ----- | | Nd 585 nm | ~585 nm | Weaker, sits on the green-red slope | | Nd 745 nm | ~745 nm | Strong and diagnostic | | Nd 810 nm | ~800-810 nm | Strong; the true minimum is near 800 nm | | Nd 870 nm | ~870 nm | At the edge of the Extended VNIR range | ### Where the reference spectrum comes from To recognize this fingerprint we need a known example to compare against. We use a laboratory spectrum of **bastnaesite** (the Nd-bearing carbonate mineral mined at Mountain Pass) from the [USGS splib07 spectral library](https://dx.doi.org/10.3133/ds1035), the standard, freely available reference library of minerals and materials measured under controlled laboratory conditions. Rather than parsing the USGS library's native format, we access it through [OpenSpecLib v0.0.6](https://github.com/null-jones/openspeclib/releases/tag/v0.0.6), an open-source project that repackages the splib07 release (a large collection of per-spectrum ASCII data files, alongside the ECOSTRESS and ECOSIS spectral libraries) as [Apache Parquet](https://parquet.apache.org/) files. That means we can pull a single reference spectrum with a couple of lines of `pandas`, with no special readers required; we simply download the pre-built Parquet file straight from the OpenSpecLib release rather than installing anything. We use the bastnaesite entry `splib07a_Bastnaesite_REE_WS320_crystl_ASDFRb_AREF`; the library also contains pure neodymium oxide and other REE minerals (monazite, samarium oxide, and more) if you want to experiment. The USGS splib07 bastnaesite reference at full 1 nm resolution (grey) and resampled to Dragonette's bands (blue); we'll cover the resampling in the "Resample a real reference spectrum" section below. Note the Nd features at the shaded windows. ## Code ### Imports and configuration ```python import io import requests import numpy as np import pandas as pd import rasterio from rasterio.warp import transform_bounds from rasterio.windows import from_bounds from scipy import ndimage from spectral.algorithms.detectors import ace as spy_ace import matplotlib.pyplot as plt # Wyvern Open Data Program STAC item (L2A surface reflectance) for the Mountain Pass scene. STAC_ITEM = ( "https://wyvern-odp.com/application/mining/" "wyvern_dragonette-002_20250422T221452_2bc40137_l2a/" "wyvern_dragonette-002_20250422T221452_2bc40137_l2a.json" ) REFLECTANCE_SCALE = 1e-4 # L2A scaled integers -> 0..1 reflectance # The ODP item is the full satellite swath, so we crop to a bounding box around the # Mountain Pass mine (min_lon, min_lat, max_lon, max_lat, WGS84) before analysis. AOI_BOUNDS = (-115.555, 35.470, -115.515, 35.492) # Pre-built Parquet asset published with the OpenSpecLib release; we download it # directly for convenience instead of installing openspeclib and converting splib07 ourselves. OPENSPECLIB_PARQUET_URL = ( "https://github.com/null-jones/openspeclib/releases/download/v0.0.6/usgs_splib07.parquet" ) BASTNAESITE_ID = "usgs_splib07:splib07a_Bastnaesite_REE_WS320_crystl_ASDFRb_AREF" ND_FEATURES = [585, 745, 810, 870] # diagnostic Nd(III) VNIR features, nm ``` After downloading the scene from the Open Data Program (see the notebook for the STAC download helper), we read the band centre wavelengths and per-band `FWHM` (Full Width at Half Maximum) straight from the GeoTIFF. Because the ODP item is the full swath, we reproject the AOI to the image CRS and crop to our AOI, then scale to reflectance: ```python def read_wavelengths_and_fwhm(dataset: rasterio.DatasetReader) -> tuple[np.ndarray, np.ndarray]: """Read band centre wavelengths and FWHM (both in nm) from a Wyvern GeoTIFF. Args: dataset: An open rasterio dataset for a Wyvern GeoTIFF. Returns: A tuple of the band centre wavelengths and the band FWHM values, both in nm. """ centres, fwhm = [], [] for band in range(1, dataset.count + 1): centres.append(float(dataset.descriptions[band - 1].split("_")[1])) fwhm.append(float(dataset.tags(band)["FWHM"])) return np.asarray(centres), np.asarray(fwhm) with rasterio.open(local_file) as src: wavelengths, fwhm = read_wavelengths_and_fwhm(src) # Reproject the AOI to the image CRS and read only that window. left, bottom, right, top = transform_bounds("EPSG:4326", src.crs, *AOI_BOUNDS) window = from_bounds(left, bottom, right, top, src.transform) window = window.round_offsets().round_lengths() cube = src.read(window=window).astype("float32") cube[cube == src.nodata] = np.nan profile = src.profile profile.update( height=cube.shape[1], width=cube.shape[2], transform=src.window_transform(window), ) cube *= REFLECTANCE_SCALE ``` ### Continuum removal An absorption feature is easiest to measure relative to the local spectral background, or "continuum". We divide each pixel's spectrum by its upper convex hull, so absorptions show up as dips below 1. This is the classic continuum-removal method of Clark & Roush (1984), [*Reflectance spectroscopy: Quantitative analysis techniques for remote sensing applications*](https://doi.org/10.1029/JB089iB07p06329) (J. Geophys. Res. 89, 6329-6340). The continuum is the convex hull fitted over the top of the spectrum, representing the background reflectance if the absorption were absent; dividing the spectrum by that hull normalizes the features to a common baseline so their depths and positions can be compared between pixels and against library spectra. The same technique underpins USGS Tetracorder and the continuum removal tools in ENVI and the EnMAP-Box. ```python def continuum_removed(wavelengths: np.ndarray, reflectance: np.ndarray) -> np.ndarray: """Continuum-remove a spectrum by dividing by its upper convex hull. Args: wavelengths: Band centre wavelengths in nm. reflectance: Reflectance values at those wavelengths. Returns: The continuum-removed spectrum (1.0 is on the continuum, lower is absorption). """ points = list(zip(wavelengths, reflectance)) hull: list[tuple[float, float]] = [] for p in points: while len(hull) >= 2: (x1, y1), (x2, y2) = hull[-2], hull[-1] if (x2 - x1) * (p[1] - y1) - (y2 - y1) * (p[0] - x1) >= 0: hull.pop() else: break hull.append(p) hull_x = [h[0] for h in hull] hull_y = [h[1] for h in hull] return reflectance / np.interp(wavelengths, hull_x, hull_y) ``` Plotting a pixel over the ore-handling area, before and after continuum removal, the Nd dips pop out at the shaded feature windows: A Mountain Pass ore pixel. After continuum removal the Nd absorptions at 585, 745 and 810 nm are clear. Notice the 810 nm dip reads a little shallow, because Dragonette's 815 nm band sits slightly off the true ~800 nm minimum. This matters less for the full-spectrum ACE detector below: it matches the shape of the whole spectrum at once, so no single band's depth decides the result. ### Resample a real reference spectrum Rather than hand-pick band positions, we compare each pixel to a real laboratory spectrum of **bastnaesite** (the Mountain Pass ore mineral) from the USGS splib07 library, packaged in [OpenSpecLib v0.0.6](https://github.com/null-jones/openspeclib/releases/tag/v0.0.6). We resample it onto Wyvern's bands by convolving with each band's response (a Gaussian approximated from the FWHM), so the reference is blurred to the same spectral resolution as the imagery. ```python def load_reference(parquet_url: str, reference_id: str) -> tuple[np.ndarray, np.ndarray]: """Load one USGS splib07 spectrum from an OpenSpecLib parquet file. Args: parquet_url: URL of the usgs_splib07.parquet OpenSpecLib release asset. reference_id: The id of the spectrum to load. Returns: A tuple of wavelengths in nm and reflectance (0..1), with bad bands set to NaN. """ raw = requests.get(parquet_url, timeout=600).content df = pd.read_parquet(io.BytesIO(raw), columns=["id", "spectral_data.values"]) values = np.asarray( df.loc[df["id"] == reference_id, "spectral_data.values"].iloc[0], dtype="float64" ) values[(values < 0) | (values > 1.5)] = np.nan # USGS bad-band flags # splib07a wavelength grid: 2151 samples, 350-2500 nm at 1 nm spacing # (measured on an ASD field spectroradiometer). wl_nm = np.linspace(0.34999999, 2.5, 2151) * 1000.0 return wl_nm, values def resample_to_bands( ref_wl: np.ndarray, ref_refl: np.ndarray, centres: np.ndarray, fwhm: np.ndarray ) -> np.ndarray: """Convolve a fine reference spectrum onto sensor bands using Gaussian responses. Args: ref_wl: Reference wavelengths in nm. ref_refl: Reference reflectance at those wavelengths. centres: Sensor band centre wavelengths in nm. fwhm: Sensor band FWHM in nm. Returns: The reference reflectance resampled to the sensor bands. Note: For simplicity this approximates each band's response as a Gaussian derived from its FWHM. A rigorous workflow would instead convolve with each band's actual relative spectral response (RSR) curve, which is published for Dragonette in the Wyvern public resources: https://github.com/Nrevyw/wyvern-public-resources/tree/main/relative-spectral-responses """ good = np.isfinite(ref_refl) wl, refl = ref_wl[good], ref_refl[good] out = np.empty(len(centres)) for i, (centre, width) in enumerate(zip(centres, fwhm)): # Gaussian approximation of the band response (see Note; real RSR curves # would be more faithful). weight = np.exp(-0.5 * ((wl - centre) / (width / 2.3548)) ** 2) out[i] = np.sum(weight * refl) / np.sum(weight) return out ref_wl, ref_refl = load_reference(OPENSPECLIB_PARQUET_URL, BASTNAESITE_ID) target = resample_to_bands(ref_wl, ref_refl, wavelengths, fwhm) ``` The USGS bastnaesite reference at full 1 nm resolution (grey) and resampled to Dragonette's bands (blue). ### Detect with ACE The [Adaptive Cosine Estimator](https://doi.org/10.1109/TSP.2005.845426) (ACE) scores how well each pixel aligns with the reference direction after the scene background has been statistically "whitened" (mean-subtracted and decorrelated using the scene covariance). ACE is insensitive to overall brightness, which makes it a robust, easy-to-threshold detector. ACE is bounded between 0 and 1 where values closer to 1 are stronger matches to the reference spectra. For a pixel spectrum $x$, target $t$, scene mean $m$ and covariance $C$, with centred vectors $\tilde{x} = x - m$ and $\tilde{t} = t - m$, the score is: $$ \mathrm{ACE}(x) = \frac{\left(\tilde{t}^{\top} C^{-1} \tilde{x}\right)^2} {\left(\tilde{t}^{\top} C^{-1} \tilde{t}\right)\left(\tilde{x}^{\top} C^{-1} \tilde{x}\right)} $$ Intuitively, $C^{-1}$ whitens the background so common, high-variance materials do not dominate, and the expression is the squared cosine of the angle between the whitened pixel and the whitened target. We use the ACE implementation in [Spectral Python (SPy)](https://www.spectralpython.net/), a common hyperspectral package, and wrap it only to drop nodata pixels and reshape the result. Because we cropped to the mine, the scene is small, so this is fast and light. ```python def detect_ace(cube: np.ndarray, target: np.ndarray) -> np.ndarray: """Adaptive Cosine Estimator detection score over a reflectance cube. The detection is delegated to Spectral Python's ace; we only drop nodata pixels and reshape the result back onto the image grid. Args: cube: Reflectance cube (bands, rows, cols) with nodata set to NaN. target: Target reflectance vector (length equal to the number of bands). Returns: A 2D array of ACE scores in 0..1, with NaN where the cube is invalid. """ bands, rows, cols = cube.shape flat = cube.reshape(bands, -1).T valid = np.isfinite(flat).all(axis=1) out = np.full(flat.shape[0], np.nan, dtype="float32") out[valid] = spy_ace(flat[valid], target) return out.reshape(rows, cols) ace = detect_ace(cube, target) ``` Mapping the ACE score next to a true-colour view, the strongest matches land squarely on the mine's ore-handling and stockpile areas: ACE detection score (right). Bright pixels are the strongest matches to the bastnaesite reference. ### Overlay the detections on the RGB preview A detection map is most useful laid over a familiar view. We threshold the ACE score, apply a little morphological cleanup to drop isolated single-pixel hits and thicken the rest for visibility, then blend the result over the true-colour composite. ```python DETECTION_THRESHOLD = 0.2 mask = np.nan_to_num(ace) > DETECTION_THRESHOLD # Additional processing: remove isolated single-pixel hits, then thicken slightly. mask = ndimage.binary_opening(mask, structure=np.ones((2, 2))) mask = ndimage.binary_dilation(mask, iterations=1) rgb = stretch(660, 550, 465) fig, ax = plt.subplots(figsize=(9, 7)) ax.imshow(rgb) overlay = np.where(mask, ace, np.nan) im = ax.imshow( overlay, cmap="turbo", vmin=DETECTION_THRESHOLD, vmax=np.nanmax(ace), alpha=0.9 ) ax.set(title=f"Nd(III) detections (ACE > {DETECTION_THRESHOLD}) over RGB") ax.axis("off") fig.colorbar(im, ax=ax, fraction=0.046, pad=0.02, label="ACE") plt.tight_layout() plt.show() ``` The 0.2 threshold is an empirical choice based on this scene's score distribution, not a magic number: background ACE scores cluster near zero (the scene median is ~0.01 and even the 99th percentile is ~0.18), so a 0.2 cutoff keeps only the extreme tail of strongest matches, which top out near 0.79. If you adapt this workflow to another scene or target, inspect the ACE histogram and tune the threshold to balance missed detections against false positives. Cleaned Nd(III) detections (ACE > 0.2) blended over the true-colour composite, concentrated on the ore-handling and stockpile areas. ### Always sanity-check A detector is only trustworthy if the pixels it flags actually look like the target. Overlaying the continuum-removed spectra of the top-scoring pixels against the reference, they share the same Nd absorptions, confirming the detection is real, not an artifact: Top-scoring pixels (orange, red = mean) reproduce the Nd absorptions of the bastnaesite reference (blue). ## Wrap-up You continuum-removed a Wyvern spectrum, resampled a real USGS bastnaesite reference to Dragonette's bands, and used ACE to map Nd-bearing material at Mountain Pass, reproducing, in miniature, the idea behind the EnMAP study. A few things to keep in mind: - **It is an exploration, not proof of REE.** Treat detections as candidates to follow up, and always sanity-check their spectra. - **Spectral resolution matters.** Dragonette's ~20-30 nm bands broaden the narrow Nd features, so measured depths are smaller than a lab or a finer-resolution sensor would see. Comparing against a full reference spectrum (as ACE does) is what keeps the detection robust. Want to go further? Try the pure neodymium-oxide reference instead of bastnaesite, threshold the ACE map to export candidate polygons, or explore other detection methods like Spectral Angle Mapper and Mixture Tuned Matched Filtering. Share what you build on our [GitHub Discussions](https://github.com/Nrevyw/wyvern-public-resources/discussions)! --- [Source: https://knowledge.wyvern.space/docs/documentation/getting_started/common_tools] # Common tools You don't need to buy anything to work with Wyvern data. Four of the six tools below are free and open source, and our tutorials cover QGIS, Python, ArcGIS Pro, and ENVI. If you're starting from scratch, pick QGIS if you'd rather not write code and Python if you would. ## QGIS Free GUI Open Source [QGIS](https://qgis.org) is free, runs on Windows, macOS, and Linux, and with the [EnMAP-Box](https://www.enmap.org/data_tools/enmapbox/) plugin it handles hyperspectral cubes without any code: spectral profiles, band math, and dimensionality reduction. Most of our [no-code tutorials](/docs/tutorials/qgis/) use it. It reads and writes nearly every raster and vector format through GDAL, and its plugin library covers most things the core application doesn't. ## ArcGIS Paid GUI Closed Source [ArcGIS](https://www.arcgis.com/home/index.html) is Esri's commercial GIS platform: ArcGIS Pro on the desktop, and ArcGIS Online and ArcGIS Enterprise for hosting and sharing. It's the standard in many government agencies and large organizations, so if your team already has licences, it's often the path of least resistance. Our [ArcGIS Pro tutorial](/docs/tutorials/arcgis_pro/) covers loading a Wyvern scene. ## ENVI Paid GUI Closed Source [ENVI](https://www.nv5geospatialsoftware.com/Products/ENVI), from NV5 Geospatial, is built specifically for spectral imagery. It has the deepest set of hyperspectral tools of anything on this page, including atmospheric correction (QUAC and FLAASH), spectral libraries, target detection, and spectral unmixing, all through a graphical interface. It's scriptable through IDL and Python, and it's a paid product. Our [ENVI tutorial](/docs/tutorials/envi/dragonette_data_in_envi) takes a Dragonette scene from import through atmospheric correction to spectral analysis. ## Python Free Code Open Source [Python](https://www.python.org/) is what most of our [tutorials](/docs/tutorials/python/) use, and it's the most flexible option on this page. A Wyvern scene loads into a NumPy array with [rasterio](https://rasterio.readthedocs.io/), band math is one line of NumPy, and [xarray](https://docs.xarray.dev/) keeps wavelengths attached to the data. The same code runs on a laptop, in a notebook on Google Colab, or in a pipeline. The [Python tutorials landing page](/docs/tutorials/python/) describes the packages we rely on. ## R Free Code Open Source [R](https://www.r-project.org/) is a language built for statistics, and it's a natural choice if your analysis ends in a model or a hypothesis test rather than a map. The current packages are [terra](https://cran.r-project.org/web/packages/terra/index.html) for rasters, [sf](https://cran.r-project.org/web/packages/sf/index.html) for vectors, and [stars](https://cran.r-project.org/web/packages/stars/index.html) for data cubes. Older tutorials often use `sp`, `rgdal`, `rgeos`, and `raster`; `rgdal` and `rgeos` were archived from CRAN in October 2023, so start new work with `terra`, `sf`, and `stars`. ## GDAL Free Code Open Source [GDAL (Geospatial Data Abstraction Library)](https://gdal.org/) is the library underneath most of the other tools here: QGIS, rasterio, and R's terra and sf all use it to read and write files. You can also use it directly. `gdalinfo` prints a scene's size, projection, and per-band wavelength tags, `gdal_translate` converts between formats and extracts bands, and `gdalwarp` reprojects and clips. It's written in C and C++, with bindings for Python and other languages. --- [Source: https://knowledge.wyvern.space/docs/documentation/getting_started/atmospheric_correction] # Atmospheric correction Before you can compare a Wyvern scene with another date, another sensor, or a library spectrum, the atmosphere has to be taken out of it. This page explains the difference between **L1B top-of-atmosphere radiance** and **L2A bottom-of-atmosphere reflectance**, which one you're likely to have, and what to do if you're starting from L1B. **L1B top-of-atmosphere (TOA) radiance** is the light signal the sensor measures, including the effects of the atmosphere such as scattering and absorption by air, aerosols, and cloud. **L2A bottom-of-atmosphere (BOA) reflectance** is a corrected product that estimates how much light the surface itself reflected, as though the atmosphere were not there. L2A is what you want for studying vegetation, water, and soil, because it's comparable between scenes, dates, and sensors. ## Which product you get Wyvern delivers both L1B and L2A. The [Open Data Program](https://opendata.wyvern.space/) publishes L2A surface reflectance unless a scene is explicitly noted otherwise, and also carries a small [L1B collection](https://wyvern-odp.com/l1b/collection.json) for anyone who wants to work from radiance. Check the `processing:level` property on the STAC item, or look for the `_l1b` or `_l2a` suffix on the item ID. The two products differ in more than the atmosphere. See [processing levels](processing_levels) for the full comparison of data type, scale factor, NoData value, and coordinate reference system. :::warning Use L2A for absorption-feature work such as spectral matching, continuum removal, and narrowband indices. In L1B radiance, the shape of the solar spectrum and atmospheric features, including the oxygen absorption near 760 nm, are larger than the surface features you are trying to measure, so those features are not recoverable without correcting the scene first. ::: ## Correcting L1B yourself If you're working from L1B, convert radiance to surface reflectance before any quantitative analysis. The first step is converting at-sensor radiance to TOA reflectance, which accounts for the solar spectrum and the sun angle. Wyvern publishes code for this in [top-of-atmosphere-processing](https://github.com/Nrevyw/wyvern-public-resources/tree/main/top-of-atmosphere-processing). Getting from there to surface reflectance means modelling what the atmosphere did to the light, which is what the tools below do. These tools read and correct Wyvern data directly: | Name | Type | Details | | ---- | ---- | ------- | | [ATCOR](https://atcor.com/) | Paid | Established atmospheric correction tool. | | [ENVI QUAC](https://www.nv5geospatialsoftware.com/docs/QUAC.html) | Paid | Built into ENVI. A scene-based approach that needs no atmospheric input parameters. | | [ACOLITE](https://github.com/acolite/acolite) | Free | Open source, developed for coastal and inland water applications. | Other tools work with Wyvern data but need manual setup to get the inputs and outputs right. These include [Py6S](https://github.com/robintw/Py6S), [SICOR](https://enmap.git-pages.gfz-potsdam.de/sicor/doc/about.html), [SIAC](https://github.com/MarcYin/SIAC), [SL2P-CCRS](https://publications.gc.ca/collections/collection_2025/rncan-nrcan/m103-3/M103-3-84-2024-eng.pdf), and [libRadtran](https://www.libradtran.org/doku.php?id=start). ## More resources NASA's Earthdata has a clear reference on [data processing levels](https://www.earthdata.nasa.gov/learn/earth-observation-data-basics/data-processing-levels). --- [Source: https://knowledge.wyvern.space/docs/documentation/agents] # Wyvern for AI Agents This page is the condensed, factual entry point for **AI agents, LLMs, and automated tools** working with Wyvern hyperspectral data. Everything here is designed to be ingested programmatically. A plain-markdown mirror of this page is served at [`/AGENTS.md`](https://knowledge.wyvern.space/AGENTS.md). ## Machine-readable endpoints | Resource | URL | | --- | --- | | Curated site map for LLMs | https://knowledge.wyvern.space/llms.txt | | Full docs text (single file) | https://knowledge.wyvern.space/llms-full.txt | | This page as plain markdown | https://knowledge.wyvern.space/AGENTS.md | | Index library (JSON) | https://raw.githubusercontent.com/Nrevyw/wyvern-public-resources/refs/heads/main/index-library/wyvern_index_library.json | | STAC catalog root (Open Data) | https://wyvern-odp.com/catalog.json | | Open Data browser (human UI) | https://opendata.wyvern.space/ | | Code, notebooks & agent skill | https://github.com/Nrevyw/wyvern-public-resources | | Sensor spectral response curves | https://github.com/Nrevyw/wyvern-public-resources/tree/main/relative-spectral-responses | ## Agent skill (recommended) A portable, tested skill teaching coding agents the full workflow — discovery, loading, band resolution, index computation, verification — lives at [`agent-skills/working-with-wyvern-data`](https://github.com/Nrevyw/wyvern-public-resources/tree/main/agent-skills) in the public resources repo. Claude Code users can copy it into `.claude/skills/`; its `SKILL.md` also works as plain context for any other framework, and its bundled STAC helper script is standard-library Python. ## Data format facts Both L2A and L1B are purchasable. Open Data publishes L2A by default; the root catalog states all imagery is L2A unless otherwise noted. A separate [L1B collection](https://wyvern-odp.com/l1b/collection.json) held 5 items as of 2026-09-19 — read the count rather than assuming it. Item IDs carry an `_l1b` or `_l2a` suffix, and `properties.processing:level` states the level directly. **L2A surface reflectance (Open Data today), license CC-BY-4.0:** - Cloud-Optimized GeoTIFF, `uint16`, **multiply by `scale` = 0.0001 to get reflectance (0–1)**; NoData = **65535** — mask *before* scaling - Projected to a **per-scene UTM zone** (`proj:epsg`, e.g. 32639); nominal GSD **5.3 m**, with per-scene `gsd` reported around 5.2 m depending on off-nadir angle - 23 bands (Standard VNIR, Dragonette-001, ~503–799 nm) or 31 bands (Extended VNIR, Dragonette-002 onward, ~445–870 nm) **L1B top-of-atmosphere radiance:** - `float32`, W·m⁻²·sr⁻¹·µm⁻¹, no scale factor; NoData = −9999 - EPSG:4326 with non-square pixels sized to 5.3 m at scene-centre latitude - Convert to ToA reflectance with [top-of-atmosphere-processing](https://github.com/Nrevyw/wyvern-public-resources/tree/main/top-of-atmosphere-processing), or [atmospherically correct](getting_started/atmospheric_correction) to surface reflectance for quantitative work. Prefer L2A for absorption-feature analysis — L1B radiance carries the solar spectrum and atmospheric features (e.g. the O₂ band near 760 nm) that swamp subtle targets. **Where band metadata lives (two authoritative per-scene sources):** ```python # 1. The GeoTIFF's own tags — already in nm, works offline with rasterio.open(path) as src: cwl = [float(src.tags(b)["wavelength"]) for b in range(1, src.count + 1)] fwhm = [float(src.tags(b)["FWHM"]) for b in range(1, src.count + 1)] nodata = src.nodata # 65535 — read it, don't hardcode ``` 2. The STAC item's COG asset `eo:bands` (`center_wavelength` / `full_width_half_max` in **µm**) and `raster:bands` (dtype / nodata / scale). ⚠️ `src.scales` is **1.0** — the GeoTIFF does *not* carry the 0.0001 reflectance scale. Take it from STAC `raster:bands`; trusting `src.scales` silently leaves data unscaled. ## STAC catalog structure Root `https://wyvern-odp.com/catalog.json` (STAC 1.0.0) has child catalogs grouping the same scenes three ways: `year/`, `application/` (agriculture, mining, coastal, forestry, …), and `product-type/{standard,extended}` — the collection files' `rel: item` links are the scenes. Item assets (keys contain spaces — quote them): | Asset key | Content | | --- | --- | | `Cloud optimized GeoTiff` | The hyperspectral imagery (with `eo:bands`, `raster:bands`) | | `Data Mask` | Valid-data mask COG | | `Pixel Quality Mask` | Quality/cloud mask COG | | `Overview image`, `Thumbnail image` | PNG previews | | `stac_metadata` | This item JSON | | `zip_file` | Everything bundled | Filterable item properties: `eo:cloud_cover`, `datetime`, `platform`, `view:sun_elevation`, `proj:epsg`, plus the item `bbox`. Note `product_type` is the literal string `"hyperspectral"` on every item and does **not** identify the band configuration — use the band count (23 = Standard, 31 = Extended) or the `product-type/{standard,extended}` catalog path. ## Code quickstart ```python import requests, rasterio import numpy as np # 1. Discover: walk product-type collections; each rel=item link is a scene. # (Send a User-Agent header: the CDN 403s Python-urllib's default.) col = requests.get("https://wyvern-odp.com/product-type/extended/collection.json").json() item_url = next(l["href"] for l in col["links"] if l["rel"] == "item") item = requests.get(item_url).json() # 2. Resolve wavelengths -> 1-based band numbers from THIS scene's metadata cog = item["assets"]["Cloud optimized GeoTiff"] bands = cog["eo:bands"] # wavelengths in µm def band_for(nm): return min(range(len(bands)), key=lambda i: abs(bands[i]["center_wavelength"] * 1000 - nm)) + 1 # 3. Load (use a windowed read for an AOI), mask, THEN scale with rasterio.open(cog["href"]) as src: red = src.read(band_for(660)).astype("float64") nir = src.read(band_for(800)).astype("float64") nodata = src.nodata # 65535, from the file scale = 0.0001 # from cog["raster:bands"]; src.scales is 1.0 red = np.where(red == nodata, np.nan, red) * scale nir = np.where(nir == nodata, np.nan, nir) * scale # 4. Compute + verify: reflectance in [0, ~1], NDVI in [-1, 1] ndvi = (nir - red) / (nir + red) ``` Index definitions with per-product-type band mappings (JSON keyed by `"Standard VNIR"` / `"Extended VNIR"`, 1-based `band_index`): [wyvern_index_library.json](https://raw.githubusercontent.com/Nrevyw/wyvern-public-resources/refs/heads/main/index-library/wyvern_index_library.json) — this is the machine-readable source of truth. The browsable [Hyperspectral Index Library](/hyperspectral_library) is generated from it. ## Recommended Python packages ```bash pip install rasterio numpy pystac requests pyproj shapely spectral matplotlib ``` | Package | Use | | --- | --- | | `rasterio` | Read/write COGs, windowed reads, reprojection — the default loader | | `numpy` | Band math, index calculation, masking | | `pystac` | Parse STAC items (`pystac.Item.from_file(url)` works against the catalog) | | `requests` | HTTP fetches (its default User-Agent avoids the CDN 403) | | `pyproj`, `shapely` | AOI geometry and CRS transforms (each scene has its own UTM zone) | | `spectral` (SPy) | Target detection (ACE, matched filter), anomaly detection (RX), SAM, MNF, unmixing | | `scikit-learn` | Classification and clustering | | `xarray` + `rioxarray`, `dask` | Labeled dimensions, time series, out-of-core processing | | `matplotlib` | Spectral plots and index maps | `pystac-client` is only useful against a STAC *API*; the Open Data catalog is static JSON, so walk `rel: item` links instead. Avoid `pysptools` — it fails to import on current Python versions; `spectral` covers the same algorithms and is maintained. ## Spectral libraries [OpenSpecLib](https://github.com/null-jones/openspeclib) (third-party, not Wyvern-maintained) amalgamates USGS Spectral Library 7, ECOSTRESS, and EcoSIS into one schema-validated structure. **Pin a release** — counts, sizes and grid layout shift between them. v0.0.6, which the REE notebook also pins, holds 32,940 spectra: 26,780 vegetation, 2,885 mineral, 1,410 water, 470 rock, 440 man-made. Download release assets directly — no install needed: ```bash BASE=https://github.com/null-jones/openspeclib/releases/download/v0.0.6 curl -sLO $BASE/usgs_splib07.parquet # 38 MB — minerals/rocks curl -sLO $BASE/wavelengths.parquet # 0.3 MB — REQUIRED for wavelengths curl -sLO $BASE/ecosis.parquet # 307 MB — optional; vegetation, large ``` Two schema details that trip up first attempts: spectra store `spectral_data.values` but **not** their wavelengths — those live in `wavelengths.parquet`, joined on `spectral_data.wavelength_grid_id`; and wavelengths are **not uniformly µm** — usgs_splib07 and ecostress grids are µm but all 43 ecosis grids (which hold the 26,780 vegetation spectra) are nm, so read `wavelength_unit` rather than hardcoding a conversion. Bad-band fill values are large negatives (e.g. `-1.23e34`) and must be masked. There is also a no-install [browser viewer](https://null-jones.github.io/openspeclib/) that can search, plot, simulate Wyvern-band downsampling, and export CSV/ENVI `.sli`. ## Spectral analysis approaches Run these on masked, scaled reflectance shaped `(rows, cols, bands)` using `spectral`: | Goal | Method | | --- | --- | | "Is material X here?" (have a reference spectrum) | **ACE** (`sp.ace`) — scale-invariant, responds to spectral shape rather than brightness | | "What's unusual here?" (no reference spectrum) | **RX** (`sp.rx`) | | Denoise / reduce bands | **MNF** (`sp.mnf`) — preferred over PCA, orders by SNR rather than variance | | Fractional abundance | endmembers (`sp.smacc`/`sp.ppi`) → `sp.unmix` | | Sharpen narrow absorptions | `sp.remove_continuum` | ⚠️ **SPy cannot consume a masked cube directly.** Wyvern swaths are rotated, so a NoData fringe is always present. `ace`, `rx`, `mnf` and `calc_stats` raise `NaNValueError`, and `spectral_angles` and `smacc` are worse — they return corrupted or all-NaN output without raising. Pass only the valid pixels as a degenerate `(N, 1, bands)` array and scatter the results back; the agent skill ships tested `valid_pixels` / `scatter_scores` helpers for exactly this. SPy has no `mtmf()` — don't report plain `matched_filter` as MTMF. For working code, thresholding, and the failure modes of each method, use the skill's [spectral-analysis reference](https://github.com/Nrevyw/wyvern-public-resources/blob/main/agent-skills/working-with-wyvern-data/references/spectral-analysis.md) rather than reimplementing from this summary. The core pipeline is **resample reference spectrum → ACE → verify hits resemble the reference**, worked end to end in the [rare earth elements notebook](https://github.com/Nrevyw/wyvern-public-resources/tree/main/tutorial-notebooks/detecting-rare-earth-elements) (neodymium at Mountain Pass) — which runs ACE on plain scaled reflectance and uses continuum removal only to inspect features in plots, not as a detection step. Read it before building a new detection workflow. **Resample reference spectra before comparing.** Library spectra are measured at 1–10 nm; Wyvern bands are 16–32 nm wide, so a reference must be convolved onto each band's response (FWHM-weighted). Naive interpolation over-weights narrow features and produces wrong scores. The agent skill ships `scripts/resample_spectra.py` for this. **Wyvern is VNIR-only, and the range depends on the product type** — Standard VNIR is 503–799 nm, Extended VNIR 445–869 nm. Check the band count before applying this table; a target at 460 nm or 860 nm is undetectable on a Standard scene: | Detectable in VNIR | Needs SWIR (not detectable) | | --- | --- | | Rare earth elements — Nd³⁺ features at 585, 745, 810, 870 nm; Standard reaches only the first two | Clays / phyllosilicates | | Ferric iron (Fe³⁺ band ~630–715 nm, both product types) | Carbonates | | Vegetation pigments, chlorophyll, stress | Hydrocarbons, most alteration minerals | | Water constituents, chlorophyll-a, turbidity | Evaporites, sulfates | **Detecting ferric iron is not the same as naming the oxide.** Hematite's discriminating minimum is ~860 nm (just inside Extended, outside Standard) and goethite's is ~900–920 nm — beyond both. Say "ferric material is present"; only Extended supports "hematite-like", and goethite cannot be located at all. Don't over-generalize to "minerals need SWIR" — REE detection in VNIR is a proven Wyvern workflow. But when a target's diagnostic features fall outside the scene's range, say so rather than presenting a weak score as a detection. **Interpreting detection results:** scores are relative to the scene, not absolute. Threshold statistically (e.g. 99.9th percentile) and state the threshold; confirm hits aren't clouds, shadows, glint, or scene edges; check that the detected pixels' mean spectrum actually resembles the target; and corroborate with a second method. "Not detected" is a valid finding. ## Common wavelength → band numbers | Use | nm | Standard VNIR | Extended VNIR | | --- | ---: | ---: | ---: | | Blue (coastal) | 445–490 | — | 1–4 | | Green | 550 | 5 | 9 | | Red | 660 | 12 | 16 | | Red edge | 712 | 17 | 21 | | Red edge | 750 | 20 | 24 | | NIR | 800 | 23 | 27 | | NIR (upper) | 870 | — | 31 | ## Pitfalls checklist 1. **Mask NoData (65535) before applying the 0.0001 scale** — scaling first turns NoData into 6.5535, which silently poisons statistics. 2. **Never hardcode band numbers across scenes** — Standard and Extended VNIR differ (band 20 is ~750 nm on Standard but ~700 nm on Extended); resolve per scene. 3. **`eo:bands` wavelengths are µm**, not nm — but the GeoTIFF's `wavelength` tags are already nm. Also: `src.scales` is 1.0, so the 0.0001 scale must come from STAC. 4. **Asset keys contain spaces** (`item["assets"]["Cloud optimized GeoTiff"]`). 5. **Each scene has its own UTM zone** — reproject to a common CRS before mosaicking or cross-scene comparison; transform lon/lat AOIs into the scene CRS before windowed reads. 6. **Python-urllib's default User-Agent gets HTTP 403** from the catalog CDN — set any custom User-Agent (the `requests` default works). 7. **The imagery host rate-limits.** Metadata/previews come from `wyvern-odp.com`, but COGs are on `wyvern-data.com`, which returns **HTTP 429** under repeated access. GDAL often reports this as `Range downloading not supported by this server!` rather than a 429 — treat it as "back off and retry" (`GDAL_HTTP_MAX_RETRY`, `GDAL_HTTP_RETRY_DELAY`), and fall back to downloading the asset if streaming keeps failing. 8. **A wavelength farther than one FWHM from the nearest band center isn't covered** by the sensor — report that rather than substituting silently. 9. **Diagnostic features outside the scene's range are undetectable** — Standard is 503–799 nm, Extended 445–869 nm. Clay, carbonate, sulfate and hydrocarbon absorptions are in SWIR; say so instead of reporting a weak match. REE and ferric iron *are* VNIR-detectable, so don't refuse those either.