# 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 - [πŸ’Ύ Loading Data](loading_data): How to load and visualize hyperspectral data in ArcGIS Pro ### 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 Tutorials on how to use Wyvern's imagery with ENVI are located at the below link: https://www.wyvern.space/wyvern-tutorial-dragonette-data-in-envi-software ## 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 - [πŸ’Ύ Loading Data](loading_data): How to load and visualize hyperspectral data in Python - [πŸ“ˆ Visualizing Spectra & Bands](visualizing_spectra): How to select pixels and visualize their spectral curve, along with 2D scatter plots. - [πŸ’Ž Detecting Rare Earth Elements](detecting_rare_earth_elements): Map neodymium (a rare earth element) with continuum removal, spectral resampling, and the ACE detector. ## 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. ### Juptyer 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. ### PySptools [PySptools](https://pysptools.sourceforge.io/index.html) is also a Python library for interacting with hyperspectral data. It contains a number of different functions ranging from machine learning to endmembers extraction algorithms. PySptools also contains a Scikit-learn interface (in alpha!), which allows it to seamlessly integrate with the [Scikit-learn](https://scikit-learn.org/stable/) ecosystem. This enables users to integrate PySptools into scikit-learn pipelines for easy ML model development. ### 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 - [πŸ’Ύ Loading Data](loading_data): How to load and visualize hyperspectral data in QGIS - [🌈 PCA Images](pca): PCA (principal component analysis) visualization of hyperspectral data with QGIS and Orfeo Toolbox - [πŸ“ˆ Visualizing Spectra](visualizing_spectra): How to select pixels and visualize their spectral curve - [πŸ”Ž Creating Indices](creating_indices): How to create hyperspectral indices in QGIS. ### 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 ## Background Multi-spectral imagery from satellites like Landsat and Sentinel-2 have typically been used to calculate 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 indices enable us to measure the presence and general overall health of vegetation within a scene. **Unfortunately, due to the width and number of bands within multi-spectral satellites, we are unable to quantify specific features of plant health like chlorophyll and carotenoid content or leaf area index.** As we can see from the images below, NDVI is still capable of reliably separating vegetation from the surrounding scene. We can also see areas of denser plant growth, and can discern lines of crops within a field. Sentinel-2 RGB Sentinel-2 NDVI Sentinel-2 RGB (true color) and calculated NDVI index We can visualize the band width of typical multi-spectral satellites with the below image. We can see that the use of a near-infrared bands allows us to effectively isolate green vegetation from soil, as vegetation has a significantly different reflectance curve. We can also see that the band widths are quite wide. For example, Sentinel-2's Band 8 (NIR, #4 in this image) has a band width of approx. 100nm. Figure 1: Wide bandwidths of multi-spectral satellites `[4]` When we compare multi-spectral bandwidths to the below image of more detailed vegetation reflectance curves, we can see that plants actually have an incredible amount of information that would otherwise be undetectable with multi-spectral satellites. It's clear that hyperspectral sensors would be able to glean significantly more information about a plant's health than multi-spectral sensors. Figure 2: Overview of crop stages & reflectance, with key areas highlighted `[5]` ## Plant Health Figure 3: Absorbance curves of plant compounds `[1]` There are many aspects of plant health that we can quantify based on the reflectance of plant foliage. - [Chlorophyll](https://en.wikipedia.org/wiki/Chlorophyll): Green pigments that allow a plant to absorb energy from light. - [Carotenoids](https://en.wikipedia.org/wiki/Carotenoid): Yellow, orange, and red pigments that have many functions within a plant. High levels of carotenoids with an absence of chlorophyll is indicative of plants that are not healthy. - [Anthocyanins](https://en.wikipedia.org/wiki/Anthocyanin): Purple, blue, or black pigments that are usually not present during the growing season. Anthocyanins are produced towards the end of summer, and are responsible (along with carotenoids) for the leaf colors we see in fall. - [Leaf Area Index (LAI)](https://en.wikipedia.org/wiki/Leaf_area_index): Dimensionless quantity that quantifies green leaf area in relation to ground surface area. - [Red Edge Position](https://en.wikipedia.org/wiki/Red_edge): Position of a rapid change in reflectance that occurs in the near-infrared range of light. Useful for crop identification, growth phase, crop monitoring, 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 process of aging in plants. Usually tied to an increase in anthocyanins (in annuals & perennials) and a decrease in chlorophyll. There are many different combinations of hyperspectral & multi-spectral bands that can be used to quantify these aspects of plant health. Below are a selection of popular indices, and their applicability to certain satellites. :::note This is just a small selection of indices! Longer lists are available from sources at the bottom of this article! ::: | Index | Applications | Sentinel‑2 | Landsat 7/8 | SuperDove | Wyvern Dragonnette | Wyvern Gen 2 | | ---- | ---- | :----: | :----: | :----: | :----: | :----: | | **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 | | | | | | **NDWI** | Sensitive to changes in canopy water content | | | | | | **CAI** | Sensitive to cellulose, useful in fire fuel conditions and crop residue monitoring | | | | | :::info Since hyperspectral imagery has narrower bandwidths, broadband indices can be replicated by combining multiple narrow bands or selecting a similar narrowband index. ::: ## Hyperspectral Indices :::note There are an incredible amount of hyperspectral indices available. With many applications including assessing plant health, urban environments, and measuring the prescence of water & snow. We will only highlight some of the most popular hyperspectral indices here, however **additional resources are listed at the bottom of this article and are worth looking into!** ::: ### RENDVI (Red Edge Normalized Difference Vegetation Index) RENDVI is the narrowband equivalent of NDVI. This index takes advantage of slight changes in sensitivity of vegetation in the red edge band (750 nm) to examine changes in foliage and senescence. RENDVI also has the advantage of not "saturating" to the extent that NDVI does (see below, some areas of plant growth in the Sentinel-2 NDVI are completely green). Sentinel-2 NDVI Wyvern RENDVI Comparison of 10m Sentinel-2 NDVI to 5m Wyvern RENDVI We are also putting the Sentinel-2 imagery at a distinct disadvantage here, as the resolution (Ground Sampling Distance) is 10m compared to Wyvern Dragonette's 5m GSD. **That being said, it's clear that RENDVI provides clear benefits to the traditional broadband NDVI**. ### VREI2 (Vogelmann Red Edge Index 2) VREI2 is a narrowband index that is effective in light or dark backgrounds [2]. It is sensitive to the combined effects of foliage cholorphyll concentration, leaf area index, and water content [3]. Wyvern RENDVI Wyvern VREI2 Comparison of 5m Wyvern RENDVI to 5m Wyvern VREI2 When comparing RENDVI to VREI2 in the above image, a number of things can be discerned. The field in the top-right is significantly greener than other fields when viewed in the VREI2 index. This field is growing maize, whereas fields in the left of the image are growing wheat and sugar beets (according to the EU 2018 crop map). Areas of more vigourous growth can also be identified in the forest at the bottom of the image in VREI2. ## Do it Yourself! Below is an example of how to calculate RENDVI from a hyperspectral image. It requires you to setup a python environment with `Numpy`, `Rasterio`, and `Matplotlib`. We highly recommend you use `Anaconda` to install these packages (`Rasterio` relies on `GDAL` which is notorious for being hard to install). Once you run this code, you should have an output that looks like the following: RGB & RENDVI Output ```python import rasterio import numpy as np import matplotlib.pyplot as plt # We're using GDAL's vsis3 driver to read from a public S3 URL. You can change this to wherever your copy # of the data is stored. INPUT_FILE_PATH = "/vsis3/wyvern-drone-data-processed/processed/SFB1516/WYVERN_SFB1516_20220726_v0_1_0_5mGSD.tiff" OUTPUT_FILE_PATH = "rendvi_output.png" # First, let's read our hyperspectral image using rasterio image = rasterio.open(INPUT_FILE_PATH) image_arr = image.read() # Indicies 23 & 19 correspond to the closest bands to the RENDVI formula. (750.179nm & 701.382nm) # Band wavelengths can be accessed by calling `image.descriptions` band_750nm = image_arr[23] band_705nm = image_arr[19] # We'll also grab a pseudo-RGB array we can plot alongside the RENDVI! # These are array indices that are close to red, green, and blue channels we would see with our eyes. rgb_array = image_arr[[12, 6, 1]] rgb_array = np.moveaxis(rgb_array, 0, -1) # Move axis just puts our RGB channel at the end of the array # Now we can use the RENDVI formula to calculate the index # 750nm - 705nm / 750nm + 705nm rendvi_array = (band_750nm - band_705nm) / (band_750nm + band_705nm) # We can now plot the array using matplotlib! fig, axs = plt.subplots(1, 2, figsize=(15, 10)) axs[0].imshow(rgb_array * 9) # Quick way to increase brightness axs[1].imshow(rendvi_array, cmap="RdYlGn") axs[0].set_title('RGB (Approx)') axs[1].set_title('RENDVI') # And now we can export as an image! plt.savefig(OUTPUT_FILE_PATH) ``` ## More Resources [Check out our deep dive into hyperspectal indices at our Index Library!](/hyperspectral_library) ## Citations `[1]` Zhou, Ning & Dong, Weiming & Mei, Xing. (2006). Realistic Simulation of Seasonal Variant Maples. 295 - 301. 10.1109/PMA.2006.29. `[2]` Vogelmann, Jim & Rock, Barrett & Moss, D.M.. (1993). Red edge spectral measurements from sugar maple leaves. Int. J. Remote Sens.. 14. 10.1080/01431169308953986. `[3]` Velichkova, Kalinka & Krezhova, Dora. (2019). COMPARATIVE ANALYSIS OF HYPERSPECTRAL VEGETATION INDICES FOR REMOTE ESTIMATION OF LEAF CHLOROPHYLL CONTENT AND PLANT STATUS. RAD Association Journal. 3. 10.21175/RadJ.2018.03.034. `[4]` 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. `[5]` [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 Before we dive head first into hyperspectral imagery, it's best to understand some of the fundamentals of geospatial processing. Wyvern's satellites produce raster imagery, which can then be used by geospatial analysts & tools like QGIS and Python to extract value from the data. ## Details In earth observation and geospatial processing, two main types of data are used: vector data and raster data. Vector data represents geographic features as points, lines, and polygons, and is stored as a series of coordinates. Vector data is often used to represent features such as roads, buildings, and land boundaries. One of the main advantages of vector data is that it can be easily scaled and manipulated, as it is defined by mathematical equations rather than pixels. Vector data is typically stored in file formats such as shapefiles or GeoJSON files. Example of vector data (grazing areas in Alberta) Raster data, on the other hand, represents geographic features as a grid of pixels, with each pixel representing a specific geographic location. Raster data is often used to represent continuous variables such as elevation, land cover, and temperature. Raster data is typically larger in size than vector data and is more difficult to manipulate, but it is well-suited for representing detailed and complex features. Raster data is typically stored in file formats such as TIFF or GeoTIFF files. Example of raster data (farmland in Germany) In earth observation, both vector and raster data are used to study a wide range of topics, including climate change, land use and land cover, natural disasters, and resource management. For example, vector data can be used to map buildings and infrastructure, while raster data can be used to study changes in land cover over time or to monitor the health of forests and other vegetation using vegetation indices like NDVI. One key difference between vector and raster data in earth observation is how they are displayed. Vector data is displayed as points and shapes, commonly colored based on vector metadata, while raster data is displayed as a grid of pixels like an image. Vector data is often used to create detailed and accurate maps, while raster data is often used to represent continuous variables such as elevation or land cover. Vector data overlaid on top of raster data In conclusion, vector and raster data are two important types of data used in earth observation and geospatial processing. Vector data represents geographic features as points, lines, and polygons, and is stored in file formats such as shapefiles or GeoJSON files. Raster data represents features as a grid of pixels (an image!) tied to a geographic location, and is stored in file formats such as TIFF or GeoTIFF files. Understanding the differences between these two types of data is essential for working with earth observation data and using GIS tools and software. ## 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 and compound reflects light differently based on it's chemical composition. Many organizations have collected, cataloged, and published these reflectance spectra for use in analysis and detection. ## USGS Spectral Library The United States Geological Service has analyzed close to 2500 individual materials for their reflectance spectra. [Their spectral library](https://crustal.usgs.gov/speclab/QueryAll07a.php) provides an easy to use tool to search & download reflectance plots from ultraviolet to short wave infrared. Screenshot of the USGS spectral library version 7 The USGS spectral library also includes reflectance spectra for a wide range of materials, including vegetation, water, and hydrocarbons. It is regularaly updated and expanded to ensure that it remains a comprehensive and valuable resource. ## 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/) that is focused on urban materials. ## ECOSTRESS Spectral Library [The ECOSTRESS Spectral Library](https://speclib.jpl.nasa.gov/) contains over 3400 natural and man made material spectra, and includes data from three other spectral libraries (Johns Hopkins, JPL, and USGS). --- [Source: https://knowledge.wyvern.space/docs/documentation/] # Knowledge Centre Welcome to Wyvern's Knowledge Centre 🌈! We want to make using hyperspectral imagery as easy and trouble-free as possible. Hyperspectral imaging has the potential to revolutionize how we look at our world, so we're as excited as you are! This area of the Knowledge Centre is focused on high-level articles around background GIS and hyperspectral knowledge, while also highlighting potential data sources, use cases, and additional resources. :::info Have questions, suggestions, 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/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/python/loading_data] # Loading & Viewing Wyvern Data in Python Output color infrared image from the Python notebook Relevant code & supporting files (anaconda environment, requirements, etc) can be found here: https://github.com/Nrevyw/wyvern-public-resources/tree/main/visualizing-wyvern-data > **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 - Ensure you have installed Anaconda/Conda - Build the environment for this notebook using the following command: ``` conda env create -f environment.yml ``` - Run this notebook in your tool of choice (Jupyter Lab, Visual Studio Code, etc), making sure you run it with the conda environment you just created! ## 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/wyvern_dragonette-001_20240930T070744_08fd7f5a.json)! ## 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 requests import rasterio import matplotlib.pyplot as plt import numpy as np ``` Then, we can define our constants. we'll use the `STAC_ITEM` object to download the STAC metadata (and subsequent URL to download the geotiff!) ```python # Feel free to modify the image that is downloaded in this section DOWNLOAD_IMAGE = True # Change to False if you have already downloaded the image! STAC_ITEM = "https://wyvern-prod-public-open-data-program.s3.ca-central-1.amazonaws.com/year/2024/wyvern_dragonette-001_20240930T070744_08fd7f5a/wyvern_dragonette-001_20240930T070744_08fd7f5a.json" LOCAL_FILE_NAME = "wyvern_dragonette-001_20240930T070744_08fd7f5a.tiff" NAN_VALUE = -9999 # Default NaN value for Wyvern L1B geotiffs ``` We also have to define a min/max normalization function. The docstring goes into good detail on why. ```python def min_max_normalize(arr: np.ndarray) -> np.ndarray: """Wyvern L1B data is Top-of-Atmosphere Radiance. These values are out of the range that matplotlib is expecting. We need to scale our image values to 0-1 or 0-255. This function does this, as well as replacing the NaN value (-9999) with an actual NaN value. Args: arr (np.ndarray): Input array requiring scaling Returns: np.ndarray: Scaled output array w/ NaN's added """ arr = np.where(arr == NAN_VALUE, np.nan, arr) # nanmin/nanmax ignore NaN values within an array while still calculating the value (max, min) return (arr - np.nanmin(arr)) / (np.nanmax(arr) - np.nanmin(arr)) ``` Now we're ready to load our STAC JSON metadata from the Open Data Program catalog! We can then look in the metadata for a URL to download the hyperspectral GeoTIFF. ```python # We're going to use the Python Requests library to download our image. Since it's a large image, # we'll use streaming/chunking to write the file directly to disk instead of holding it in memory. # Feel free to not run this section if you have already downloaded the data. print("Loading STAC Item!") stac_item_response = requests.get(STAC_ITEM) stac_item_response.raise_for_status() # Will raise an error if we have any issues getting the STAC item stac_item = stac_item_response.json() print(f"Successfully loaded STAC Item!\nSTAC Item ID: {stac_item['id']}") download_url = stac_item["assets"]["Cloud optimized GeoTiff"]["href"] local_file_name = download_url.split("/")[-1] print( f"Downloading GeoTIFF from STAC Metadata!\nDownload url: {download_url}]\n" f"Downloading to local file named: {local_file_name}" ) if DOWNLOAD_IMAGE: with requests.get(download_url, stream=True) as r: r.raise_for_status() with open(local_file_name, "wb") as f: for chunk in r.iter_content(chunk_size=8192): f.write(chunk) print("Download completed!") ``` Once we have the image downloaded, we can load it using [rasterio](https://rasterio.readthedocs.io/en/stable/api/index.html). ```python # Now, let's load in our downloaded image & print some metadata about the image. image_file = rasterio.open(LOCAL_FILE_NAME) image_arr = image_file.read() print(f"Image shape: {image_arr.shape}") print(f"Image bands: {', '.join(image_file.descriptions)}") ``` Now that we have the image array loaded, we're good to plot it using [matplotlib](https://matplotlib.org/)! ```python # First, we select our R, G, and B channels from our hyperspectral image array, and then we shift the array shape # from [BAND, X, Y] to [X, Y, Band] which matplotlib is expected rgb_arr = image_arr[[10, 4, 0], :, :].swapaxes(0, -1) # Then we can normalize our image and plot the result! # 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(min_max_normalize(rgb_arr)) # For fun, let's also export this as an image! plt.tight_layout() plt.savefig("wyvern_rgb.png") ``` Output true-color image from the Python notebook We can also render the image in color infrared by selected a different set of bands from the image. ```python # Now, let's do the same with color infrared! cir_arr = image_arr[[19, 10, 4], :, :].swapaxes(0, -1) # Then we can normalize our image and plot the result! fig, ax = plt.subplots(1, 1, figsize=(15, 10)) ax.imshow(min_max_normalize(cir_arr)) # For fun, let's also export this as an image! plt.tight_layout() plt.savefig("wyvern_cir.png") ``` 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 image easy to use in your software of choice. This section of the knowledge centre is focused on easy to use and follow tutorials/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/datasets_and_use_cases/monitoring_carbon] # Monitoring Carbon Capture Carbon credits are a crucial tool used to reduce greenhouse gas emissions contributing to climate change. To validate carbon credits and manage carbon at scale, satellites will need to play a larger role because they can provide objective information over huge areas at an affordable price. Hyperspectral data is well suited to more accurate and efficient validation, and could be the key to expanding the carbon credit market worldwide. Hyperspectral imagery can sometimes infer net change in carbon in certain ecosystems. It can also directly detect greenhouse gas emissions. Here, we describe how the technology could be used by analytics companies to improve carbon credit accountability. ## Background Carbon credits are important because they provide a financial incentive for businesses and individuals to reduce their greenhouse gas emissions. By creating a market for carbon emissions, carbon credits give emitters a way to offset their emissions by funding projects that reduce carbon emissions elsewhere. This can help to reduce the overall level of carbon dioxide in the atmosphere and mitigate the effects of climate change. In addition, the use of carbon credits allows for the implementation of a cap and trade system, where a government or regulatory body sets a limit on the total amount of carbon that can be emitted. Emitters are then issued a certain number of carbon credits, which they can trade with each other if they need to exceed their limit. This creates a financial incentive for emitters to reduce their own emissions, as they can sell any excess credits they have for a profit. Many countries have instituted some form of cap and trade system. Carbon offsets can be a useful tool for individuals and businesses looking to reduce their own carbon footprint, but it is important to ensure that the offset projects are legitimate and effective. This can involve verifying the veracity of the claims made by the project developers, and ensuring that the offsets are additional, meaning that the emissions reductions would not have occurred without the funding provided by the offset purchase. ### Agricultural Land Management Agricultural Land Management (ALM) centers around sequestering carbon in the soil and improving farming practices. Some of the most common practices within ALM are as follows: - Improving planting & harvesting practices through cover crops & rotation - Reducing soil tillage - Improving irrigation practices - Improving grazing practices - Reducing fertilizer usage :::note Most of the approved methods of assigning carbon credits in agriculture centers around *validating* practices and *not* direct measurement. Direct measurement of carbon using imaging is still not proven. However, validating agricultural practices is a proven correlation to carbon sequestration potential, and hyperspectral imaging will greatly enhance the effectiveness. ::: Verra, a standards setting organization, provides an extremely detailed methodology for improved ALM [here](https://verra.org/methodologies/vm0042-improved-agricultural-land-management-v2-1/), which details measurement methods & formulae for calculating carbon credits. ### Improved Forest Management Improved Forest Management (IFM) is a number of different approaches that focus on conservation and protection of forests to increase the amount of carbon sequestered over a forest's lifespan. Common practices include: - Extending rotation age of forests prior to harvesting - Fire Management - Rehabilitation of previously logged forests :::note Similar to ALM, direct measurement of carbon is not yet possible with imaging. However there is strong relationships between forestry practices, vegetation condition, and practice validation that can be directly correlated to carbon credits and sequestration potential. ::: Verra also has methodologies for improved forest management [here](https://verra.org/methodologies-main/). ## Wyvern's Solution Hyperspectral imaging is uniquely suited to validating practice management for ALM and IFM. Because hyperspectral imaging provides more spectral information, additional detail can be extracted from an image to classify a characteristic of interest. ### Soil Tillage A key aspect of ALM carbon credits is the reduction of soil tillage. This can easily be monitored using hyperspectral imagery to detect the presence of lignin and cellulose, and the abscence of healthy plant compounds like chlorophyll and carotenoids. Tilled soil will contain less lignin and cellulose as compared to untilled soil. One such measure is the Lignin Cellulose Absorption Index, which is sensitive to the absorption of lignin and cellulose in the SWIR spectrum `[1]`. When we apply the **LCAI** equation to hyperspectral imagery, areas that contain large amounts of dead plant matter (and therefore lots of cellulose & lignin without chlorophyll) will show up as positive values. Healthy vegetation and dirt will present as negative values. This can be observed in the below image, where fields containing large amounts of dead plant material show up as bright green, and healthy fields & dirt present as red & yellow. Wyvern RGB Approx. Wyvern LCAI Comparison of pseudo-RGB imagery to LCAI for measuring crop residue ### Cover Crops & Plant Growth Detection of cover crops on fields can easily be accomplished by using hyperspectral indices. Additionally, hyperspectral imagery can be used to easily train classification algorithms to detect what types of cover crops are growing on a field. A number of hyperspectral indices are well suited for assessing cover crops & crop health. These include RENDVI, and VREI2. RENDVI is useful to detect the presence of healthy vegetation, as non-vegetation is easily filtered out. Wyvern RGB Wyvern RENDVI Comparison of 5m Wyvern RGB to 5m Wyvern RENDVI As we can see above, RENDVI is effective at isolating plant growth from fallow fields and infrastructure. Incorporating the use of RENDVI with a time-series analysis of fields over winter would provide an excellent solution to tracking cover crop usage across large areas. ## More Resources [Verra Methodology for Improved Agricultural Land Management](https://verra.org/methodologies/vm0042-improved-agricultural-land-management-v2-1/) ## Citations `[1]` Daughtry, Craig & Hunt, E.R. & Doraiswamy, Chandraprabha & McMurtrey, J.E.. (2005). Remote Sensing the Spatial Distribution of Crop Residues. Agronomy journal. 97. 10.2134/agronj2003.0291. --- [Source: https://knowledge.wyvern.space/docs/documentation/getting_started/hyperspectral_basics] # Hyperspectral Basics Hyperspectral imaging is a relatively new technology to be commercialized. It is a type of imaging that has a wide range of potential use cases, from food processing to observing distant stars. We will be looking at hyperspectral imaging from an earth observation (EO) lens, as we operate hyperspectral satellites! It's best to start with the electromagnetic spectrum (what we're sensing when we talk about imagery), and then different types of earth observation imaging of the electromagnetic spectrum, and you'll see how hyperspectral fits in! ## Electromagnetic Spectrum The electromagnetic spectrum refers to the full range of electromagnetic waves that are emitted by various sources. These waves have different wavelengths and frequencies, and they can be used for a variety of purposes. The electromagnetic spectrum is often divided into different categories, such as visible and near-infrared light. Visible light is the part of the electromagnetic spectrum that we can see with our eyes, while near-infrared light is just beyond the range of human vision. Both of these types of light can be used for remote sensing, as they can be detected by specialized instruments and used to gather information about the Earth and its environment. The electromagnic spectrum, with a zoomed in view of visible light\'s spectrum. `[4]` Wikipedia provides a considerably [more detailed explanation](https://en.wikipedia.org/wiki/Electromagnetic_spectrum) about the electromagnetic spectrum than we can provide here. In short, satellites carry instruments to sense light reflected by objects on earth, which can then be used to infer information about an objects properties. ## Imagery Types Now that we know that we're sensing the reflectance of electromagnetic radiation from objects on earth, we can dig into ways we image this electromagnetic radiation. While there are many different approaches, we'll be focusing on three common types of imaging in earth observation (panchromatic, multi-spectral, and hyperspectral). These types of imaging are differentiated based on how they sense the electromagnetic spectrum. A visual example of Panchromatic, Multi-Spectral, and Hyperspectral (with dimensionality reduction) imagery; from left-to-right. ### Panchromatic Panchromatic imagery is typically used to provide a high-resolution image, which is then combined with other, low resolution bands in a process known as "pansharpening". Generating a high-resolution image typically requires **more** light than a low-resolution image, so instead of filtering out light into "red" or "blue" bands, all light is passed directly to the sensor. The downside of a panchromatic image is that it contains no wavelength specific ("color") information, instead only showing the relative brightness of objects. Think of it like a black and white image! Greyscale example of a panchromatic band Since panchromatic bands are typically designed to pansharpen lower resolution wavelength-specific bands from multi-spectral or hyper-spectral imagers, a panchromatic band will collect light from a similar range of wavelengths that the instrument is sensitive to. :::danger Pansharpening bands outside of wavelengths the panchromatic imager is sensitive to carries the risk of inaccurate results! ::: ### Multi-spectral Multi-spectral imagery differs from panchromatic imagery in the sense that it is collecting **wavelength-specific** information. It uses fairly wide bandwidths (think of it like collecting most of the light we would consider "blue" or "red"). Typically MSI collects information about red, green, blue, and near-infrared wavelengths; with some imagers also collect short-wave infrared (SWIR), and the area between red and near-infrared ("red edge"). Since multi-spectral imagers typically use wide bandwidths, there is the potential for "bleed-over", where light from one band is also picked up by another band. Figure 1, below, shows a good example of wide band-widths and "bleed-over". The "Blue" and "Green" bands can be seen to have significant overlap of their sensitivity. Figure 1: Spectral response and band overlap ("bleed-over") of a typical, 4 band multispectral satellite. Below in Figure 2 is an example of 3 commmonly imaged items; soil, vegetation (crops, forest, etc), and water. Multi-spectral bandwidths are represented by the grey bars. It can be seen that typical multi-spectral bandwidths do not cover the entire spectrum, and lack considerable information. That being said, multi-spectral bands are still capable of differentiating between commmonly imaged spectra, specficially due to the considerably different reflectance of vegetation in the near-infrared and mid-infrared ranges. Figure 2: Commonly imaged spectra `[2]` ### Hyperspectral Hyperspectral imagery goes a step further than multi-spectral imagery. Instead of a small number of bands with wide bandwidths, hyperspectral imagery has a large number of bands with small band widths. This allows a considerably more detailed analysis of reflectance spectra. Figure 3 details a good example of many, small band widths covering the visible to near-infrared (VNIR) spectrum. Figure 3: Example spacing of hyperspectral bands across the VNIR spectrum There are significant advantages to hyperspectral bands for measuring everything from vegetation health, soil quality analysis, to greenhouse gas (GHG) detection. Wyvern's first generation of satellites will provide bands in the VNIR spectrum, while our second and third generation satellites will bring detailed hyperspectral monitoring into the SWIR spectrum. There are also disadvantages to hyperspectral imagery. Due to the increased number of bands, hyperspectral imagery has considerably bigger file sizes. This can result in data storage challenges, as well as increased processing time for analysis of imagery. Figure 4: Example of hyperspectral reflectance curves for various vegetation `[3]` In Figure 4, we can see that various types of vegetation have different reflectance curves. With hyperspectral band widths, we can easily start to differentiate vegetation based on how it reflects light. With multi-spectral imagery's wide band widths, we may lose some of the nuance of the data, making differentiation of different vegatation types more difficult. This also applies to other substances like soil, minerals, water, and more! ## 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 ## AVIRIS / AVIRIS-NG The Airborne Visible / Infrared Imaging Spectrometer is a remote sensing instrument that is mounted on an aircraft and used to collect spectral data of the Earth's surface. This data can be used to study the composition and characteristics of the Earth's surface, including vegetation, minerals, and water. AVIRIS uses a range of wavelengths in the visible and infrared spectrum to provide detailed information about the materials it measures. This information can be used in a variety of applications, such as mapping and monitoring the environment, studying the impacts of natural disasters, and managing natural resources. It is a project by NASA's Jet Propulsion Laboratory, and has been operating since 1986 across the globe. AVIRIS has two data portals to easily access datasets: [AVIRIS Data Portal (2006-2021)](https://aviris.jpl.nasa.gov/dataportal/): Contains imagery from the first AVIRIS sensor. This data portal also contains significantly more data since it comprises roughly 15 years of flights, however it is mostly localized to Canada & USA imagery. [AVIRIS-NG Data Portal (2014-2021)](https://avirisng.jpl.nasa.gov/dataportal/): Contains imagery from the AVIRIS-NG sensor. AVIRIS-NG overflights include Europe, India, and the arctic. ## EnMAP The [Environmental Mapping and Analysis Program](https://www.enmap.org/) is a hyperspectral imagery mission by Germany which delivers imagery at a 30m GSD, with 230 bands from 420nm to 2450nm (VNIR & SWIR). It was launched in April 2022, and started routine operations in November 2022. EnMAP requires registration to access data. ## NEON The [National Ecological Observatory Network](https://www.neonscience.org/) is operated by the National Science Foundation, and is designed to collect ecological data across the United States. The NEON program collects many different types of data, including [high resolution hyperspectral imagery](https://www.neonscience.org/data-collection/airborne-remote-sensing) of the VNIR and SWIR spectrums. --- [Source: https://knowledge.wyvern.space/docs/tutorials/python/visualizing_spectra] # Visualizing Spectra & Plotting One of the outputs from this tutorial, a pixels spectral curve! :::info Relevant code & supporting files (anaconda environment, requirements, etc) can be found here: https://github.com/Nrevyw/wyvern-public-resources/tree/main/visualizing-wyvern-data ::: 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 - Ensure you have installed Anaconda/Conda - Build the environment for this notebook using the following command: ``` conda env create -f environment.yml ``` - Run this notebook in your tool of choice (Jupyter Lab, Visual Studio Code, etc), making sure you run it with the conda environment you just created! ## 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/wyvern_dragonette-001_20240930T070744_08fd7f5a.json)! ### Spectral Plots These plots show an individual pixel's reflectance/radiance 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 showcases 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 requests import rasterio import matplotlib.pyplot as plt import numpy as np from matplotlib.patches import Rectangle ``` Then, we can define our constants. we'll use the `STAC_ITEM` object to download the STAC metadata (and subsequent URL to download the geotiff!) ```python # Feel free to modify the STAC item/image that is downloaded in this section DOWNLOAD_IMAGE = True STAC_ITEM = "https://wyvern-prod-public-open-data-program.s3.ca-central-1.amazonaws.com/year/2024/wyvern_dragonette-001_20240930T070744_08fd7f5a/wyvern_dragonette-001_20240930T070744_08fd7f5a.json" LOCAL_FILE_NAME = "wyvern_dragonette-001_20240930T070744_08fd7f5a.tiff" ``` Now we're ready to load our STAC JSON metadata from the Open Data Program catalog! We can then look in the metadata for a URL to download the hyperspectral GeoTIFF. ```python print("Loading STAC Item!") stac_item_response = requests.get(STAC_ITEM) stac_item_response.raise_for_status() # Will raise an error if we have any issues getting the STAC item stac_item = stac_item_response.json() print(f"Successfully loaded STAC Item!\nSTAC Item ID: {stac_item['id']}") download_url = stac_item["assets"]["Cloud optimized GeoTiff"]["href"] local_file_name = download_url.split("/")[-1] print( f"Downloading GeoTIFF from STAC Metadata!\nDownload url: {download_url}]\n" f"Downloading to local file named: {local_file_name}" ) if DOWNLOAD_IMAGE: with requests.get(download_url, stream=True) as r: r.raise_for_status() with open(local_file_name, "wb") as f: for chunk in r.iter_content(chunk_size=8192): f.write(chunk) print("Download completed!") ``` Once we have the image downloaded, we can load it using [rasterio](https://rasterio.readthedocs.io/en/stable/api/index.html). ```python # Now, let's load in our downloaded image & print some metadata about the image. image_file = rasterio.open(LOCAL_FILE_NAME) image_arr = image_file.read() # Replace -9999 w/ NaNs image_arr = np.where(image_arr == image_file.nodata, np.nan, image_arr) print(f"Image shape: {image_arr.shape}") print(f"Image bands: {', '.join(image_file.descriptions)}") ``` # Spectral Curve Plot We'll need our min/max normalization function here too so we can render an RGB image. ```python def min_max_normalize(arr: np.ndarray) -> np.ndarray: """Wyvern L1B data is Top-of-Atmosphere Radiance. These values are out of the range that matplotlib is expecting. We need to scale our image values to 0-1 or 0-255. Args: arr (np.ndarray): Input array requiring scaling Returns: np.ndarray: Scaled output array """ # nanmin/nanmax ignore NaN values within an array while still calculating the value (max, min) return (arr - np.nanmin(arr)) / (np.nanmax(arr) - np.nanmin(arr)) ``` Now, we're ready to plot our spectral curve. We'll select an arbritrary X/Y 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_X, PIXEL_Y = 1600, 1900 # Point on image to plot, feel free to change! pixel_spectra = image_arr[:, PIXEL_X, PIXEL_Y] # Grabbing pixel from image array band_wavelengths = [ float(x.split("Band_")[1]) for x in image_file.descriptions ] # Calculating wavelengths for line plot # 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("Top of Atmosphere Radiance") # Secondly, plotting a buffered area around our pixel to visualize what we're looking at! axs[1].imshow( min_max_normalize( image_arr[ [10, 4, 0], # R, G, B PIXEL_X - 200 : PIXEL_X + 200, # Buffering around point PIXEL_Y - 200 : PIXEL_Y + 200, # Buffering around point ].swapaxes( 0, -1 ) # Moving [Band, X, Y] to [X, Y, Band] shape ) ) # 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 Note the large downward spike near 760nm. This is an area of the electromagnetic spectrum where oxygen in the earth's atmosphere absorbs light. This means that our satellite receives less information on this band. This is why it is **very** important to perform atmospheric corrections prior to doing any analysis on Wyvern L1B imagery! # 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 X_IDX = 20 Y_IDX = 10 # Ravel turns our 2D X-Y image array into a single dimensional array band_x = image_arr[X_IDX].ravel() band_y = image_arr[Y_IDX].ravel() # Here we're filtering out our NaNs from the array band_x = band_x[~np.isnan(band_x)] band_y = band_y[~np.isnan(band_y)] # 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. heatmap, xedges, yedges = np.histogram2d( band_x, band_y, bins=200, range=[[40, 150], [40, 150]] ) 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(image_file.descriptions[X_IDX]) plt.ylabel(image_file.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! pixel_rect = Rectangle( (pixel_spectra[X_IDX] - 5, pixel_spectra[Y_IDX] - 5), 10, 10, 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! X_IDX = 22 Y_IDX = 0 # Ravel turns our 2D X-Y image array into a single dimensional array band_x = image_arr[X_IDX].ravel() band_y = image_arr[Y_IDX].ravel() # Here we're filtering out our NaNs from the array band_x = band_x[~np.isnan(band_x)] band_y = band_y[~np.isnan(band_y)] # 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. heatmap, xedges, yedges = np.histogram2d( band_x, band_y, bins=200, range=[[40, 150], [40, 150]] ) 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(image_file.descriptions[X_IDX]) plt.ylabel(image_file.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! pixel_rect = Rectangle( (pixel_spectra[X_IDX] - 5, pixel_spectra[Y_IDX] - 5), 10, 10, 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. We can see a large group of high-reflectance pixels in the centre-ish portion of each plot, with values at a similar level for each band. This group is the dirt/sand around the irrigated fields. Our X-axis on both images is in the near-infrared range, and we can see groups of pixels that have a high X (near-infrared) reflectance, but a low blue/green/red reflectance. These pixels are vegetation, since vegetation has high reflectance in the near-infrared portion of the spectrum, and low reflectance in the blue/green section of the spectrum! We can see some distinct groups in this area (feel free to zoom the plot in to this area!), these groups are likely different species of crops, or crops in different growth stages. How cool! 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/#/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/datasets_and_use_cases/moisture] # Soil & Vegetation Moisture Measuring soil and vegetation moisture is an important aspect of modern agriculture, as it can affect crop growth and yield, and can also be used to optimize irrigation and fertilization practices. There are several methods that are commonly used to measure soil and vegetation moisture in agriculture, each with its own advantages and disadvantages. High resolution hyperspectral imaging in the SWIR spectrum promises to add another effective tool to the tool box for assessing moisture. Wyvern's second generation of satellites will be able to deliver high resolution & revisit rates, enabling farmers to effectively track soil moisture during early season planning & seeding, and vegetation moisture throughout the rest of the growing season. ## Background There are two main types of measurements for soil moisture, **direct** (field sampling), and **indirect** (sensors, remote sensing, satellites, aerial imagery). Direct measurement methods involve physically extracting a sample of soil and measuring its moisture content directly, while indirect measurement methods use other properties of the soil that are related to soil moisture content, such as soil tension, soil reflectance, or soil temperature, to infer the moisture content. Direct measurement methods are generally considered to be more accurate than indirect measurement methods, but they can be more time-consuming and expensive to perform. Indirect measurement methods can be used to quickly and cheaply measure soil moisture in large areas of land, but they may not be as accurate as direct measurement methods. In the case of optical, indirect measurement methods, we use wavelengths from the visible to shortwave infrared parts of the electromagnetic spectrum (700nm to 2200nm). Certain wavelengths in this range will change in response to the moisture content of soils and leaves, which we can then use to indirectly estimate moisture content. :::info Satellites are unable to measure soil moisture once vegetation is covering the soil. Once this happens, we can still assess plant health by measuring moisture within a plant's leaves! ::: ### Current Remote Sensing Methods Remote sensing satellites are widely used to measure soil moisture, as they can cover large areas of land quickly and provide information on the distribution of moisture across an entire field or region. There are several remote sensing satellites that are currently in operation that are specifically designed for soil moisture measurements, including SMOS, SMAP, MODIS, and Sentinel-1. [SMOS (Soil Moisture and Ocean Salinity)](https://www.eoportal.org/satellite-missions/smos#smos-soil-moisture-and-ocean-salinity-mission) is a European Space Agency (ESA) satellite that was launched in 2009. SMOS is equipped with a L-band radiometer that measures the brightness temperature of the surface, which can be used to infer the soil moisture content. SMOS has a spatial resolution of about 35 km, which means that each pixel in the image represents an area of about 35 km x 35 km on the ground. [SMAP (Soil Moisture Active Passive)](https://www.jpl.nasa.gov/news/new-satellite-data-will-help-farmers-facing-drought/) is a satellite developed by NASA, launched in 2015. It uses a L-band radar and a L-band radiometer to measure soil moisture at a resolution of around 3 km, however its data is heavily affected by the vegetation cover. [MODIS (Moderate Resolution Imaging Spectroradiometer)](https://modis.gsfc.nasa.gov/about/media/modis_brochure.pdf) is a sensor aboard NASA's Terra and Aqua satellites that have been launched in 1999 and 2002 respectively, it provides data with moderate spatial resolution of around 250m and high temporal resolution with daily coverage. MODIS uses infrared and visible wavelengths to estimate surface temperature, which can be used to infer soil moisture content. [Sentinel-1](https://sentinels.copernicus.eu/copernicus/sentinel-1) is a radar imaging satellite that is part of the European Space Agency's Copernicus program. It provides data with high temporal resolution with a revisit time of 6-12 days, but its spatial resolution is quite low, around 20m-100m. Sentinel-1 can provide data even in cloudy conditions, and the radar signal can penetrate through vegetation, making it useful for monitoring soil moisture in vegetated areas. Although these remote sensing satellites can provide valuable information on soil moisture, they have some downsides. One of the main downsides is that they have low ground sample distance (GSD), which means that each pixel in the image represents a relatively large area on the ground. For example, SMOS has a GSD of about 35 km, which means that each pixel in the image represents an area of about 35 km x 35 km on the ground. This can make it difficult to obtain accurate soil moisture measurements for small areas, such as individual fields or small plots. Another downside is that the revisit rate is low, which means that a given location is only covered by the satellite every several days. This can make it difficult to observe changes in soil moisture in real-time or over short timescales. ### How Wyvern is addressing these challenges Wyvern will be launching a fleet of innovative satellites that utilize a deployable design by 2030. This allows us to reduce launch costs while still providing high-resolution hyperspectral imaging capabilities in the VNIR and SWIR regions. With the ability to measure various parameters such as vegetation health, water content and stress levels, as well as monitor soil moisture at a high spatial resolution, our technology can provide valuable insights for customers in industries such as agriculture, forestry, and land management. These insights can help optimize crop yields, reduce water usage, and improve overall efficiency in land management. ## Measuring Moisture with SWIR Short wave infrared (SWIR) imagery can be used to measure soil and vegetation moisture by analyzing reflectance in the SWIR wavelength range (900nm to 1700nm). The theory behind this method is based on the fact that the absorption of water within vegetation and soil results in changes to reflectance. ### Soil Moisture Water in the soil absorbs energy in the SWIR wavelength range, which causes a decrease in the reflectance of the soil. This decrease in reflectance can be used to infer the soil moisture content. The amount of the reflectance change is related to the amount of water in the soil and therefore, the soil moisture content can be estimated by measuring the reflectance of the soil in the SWIR wavelength range. Reflectance of a soil sample with varying degrees of moisture content `[1]` The relationship between reflectance and soil moisture content is not linear and can be affected by other soil properties such as soil type, soil texture, soil organic matter and others. Therefore, it is essential to have a good understanding of the soil characteristics of the area being measured and to calibrate the SWIR data with in-situ measurements of soil moisture to obtain accurate soil moisture estimates. ### Vegetation Moisture Vegetation is composed of cellulose, lignin, and other organic compounds that absorb energy in the SWIR wavelength range, the amount of absorption of energy is also related to the amount of water in the vegetation, meaning that the more water present the more energy gets absorbed (and therefore less light is reflected). By measuring the reflectance of vegetation in the SWIR wavelength range, it is possible to infer the vegetation moisture content. Increased canopy water content (Cw) results in decreased reflectance at 1240nm for multiple leaf area indices (simulated) `[2]` However, similar to soil moisture measurements, the relationship between reflectance and vegetation moisture content is not linear and can be affected by other vegetation properties such as species, age, leaf angle, canopy structure, etc. Therefore, it is essential to have a good understanding of the vegetation characteristics of the area being measured, and to calibrate the SWIR data with in-situ measurements of vegetation moisture to obtain accurate vegetation moisture estimates. ### Moisture Indices We can harness the change in reflectance due to water content by isolating parts of the VNIR and SWIR light spectrum that are most sensitive to this change, and developing equations that isolate this change in reflectance from other factors (chlorophyll, different materials, etc). A number of indices have been developed for measuring the water content in soils and vegetation, and include Normalized Difference Water Index (NDWI), Normalized Multi-band Drought Index (NMDI), and Normalized Difference Infrared Index (NDII). More information around index formulas and supporting details is available on our [🌟 Index Library](/hyperspectral_library). One key index, is the **Normalized Multi-band Drought Index** developed by Lingli Wang and John J. Qu. It was originally developed for use with MODIS bands, however can be applied to any hyperspectral imagery with similar band centers and band widths. Analysis of the NMDI index showed that it is effective at estimating soil and vegetation moisture in **bare soil, weekly vegetated areas, or heavily vegetated areas with a leaf area index of at least 2 [3].** Areas with moderate vegetation cause the signal from the vegetation and soil to become mixed, resulting in inaccurate results. NMDI is designed to be a soil drought measure, or vegetation drought measure; it is not able to be both at the same time. NMDI is also unable to **directly** measure moisture content of soils accurately. Below is an example of NMDI applied to a hyperspectral image. Vegetated areas present as blue, whereas bare soil and soil with light vegetation presents as brown. Areas of high and low moisture can be discerned within farmers fields. Wyvern RGB (Approx.) Wyvern NMDI Comparison of RGB imagery to NMDI. Areas of less growth in the centre field can be seen to have similar moisture levels to the rest of the field. ## More Resources [Analysis of Soil Moisture Indices](https://www.researchgate.net/publication/356348833_Optical_Remote_Sensing_Indexes_of_Soil_Moisture_Evaluation_and_Improvement_Based_on_Aircraft_Experiment_Observations) ## Citations `[1]` Philpot, William. (2010). Spectral Reflectance of Wetted Soils. 10.13140/2.1.2306.0169. `[2]` Lingli Wang, John J. Qu, Xianjun Hao & Qingping Zhu (2008) Sensitivity studies of the moisture effects on MODIS SWIR reflectance and vegetation water indices, International Journal of Remote Sensing, 29:24, 7065-7075, DOI: 10.1080/01431160802226034 `[3]` Wang, Lingli & Qu, John. (2007). NMDI: A normalized multi-band drought index for monitoring soil and vegetation moisture with satellite remote sensing. Geophysical Research Letters - GEOPHYS RES LETT. 34. 10.1029/2007GL031021. --- [Source: https://knowledge.wyvern.space/docs/documentation/getting_started/hyperspectral_image_formats] # Hyperspectral Image Formats Hyperspectral imagery is similar to most other forms of geospatial imagery, commonly called [rasters](https://desktop.arcgis.com/en/arcmap/10.3/manage-data/raster-and-images/what-is-raster-data.htm). The main difference between images generated by satellites like Landsat & Sentinel, and hyperspectral imagery is the **number** of bands within an image. The below table shows common sources of imagery and their bands. You can see our satellites will have a **ton** more bands than common multi-spectral satellites. Take note, this also means larger file sizes! | Satellite | Satellite Type | # of Bands | Common Format | | --------- | -------------- | ---------- | ------------- | | Sentinel-2 | Multi-spectral | 13 | JPEG2000, GeoTIFF | | Landsat | Multi-spectal | 11 | GeoTIFF (`.tiff`) | | Sentinel-1 | SAR | ~4 | GeoTIFF (`.tiff`) | | **Wyvern Dragonnette**| **Hyperspectral** | **23 - 31** | **GeoTIFF (`.tiff`)** | ## Background Raster imagery can be visualized as a data "cube", with 3 axis: `X (longitude)`, `Y (latitude)`, and `Z (bands)`. When we take a "slice" across the cube by selecting a single band, we get a greyscale image. Taking three "slices" and assigning each slice to a color channel (red, green, blue) would net us a color image! Three dimensional hyperspectral data "cube" `[1]` Regardless of file format, a raster is a raster. Once loaded into your software of choice, every raster is a multi-dimensional array of pixel values. When in doubt, try to load it using GDAL. Chances are it will be able to read the file. ## Common Formats We'll only cover the most common formats, however there is truly a **ludicrous** amount of image formats out there. For more details on formats (and a set of tools that can do it all!) check out GDAL's documentation on [raster drivers](https://gdal.org/drivers/raster/index.html). ### GeoTIFF `.tiff` GeoTIFF is a public domain standard developed by Dr. Niles Ritter at NASA JPL. It is extremely popular, and has many different options and variants designed to optimize use of rasters in different settings (BigTIFF, Cloud Optimized GeoTIFF). A single GeoTIFF file contains all the necessary geospatial metadata to link the raster image data to it's location in the world. [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 & .hdr` ENVI is a file format developed by NV5 for their software `ENVI`, and is popular with hyperspectral files. It consists of two parts: the actual image file, and an accompaning ASCII text header file (usually with the `.hdr` extension). GDAL (and thus rasterio and many other tools) is able to read and write ENVI files. [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 popular format in scientific computing and is maintained by the non-profit HDF Group. HDF5 files are structured in a filesystem-like format, where named groups can contain datasets, or more groups. Metadata is attached to individual datasets and groups. HDF5 files can be read using a wide variety of packages and programming languages. NetCDF4 is also based off of HDF5, and deserves a special mention since it can be processed using the powerful `xarray` library in python. [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 [Ultimate 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/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 An Earth observation sensor reads the intensity of the [electromagnetic spectrum](hyperspectral_basics?id=electromagnetic-spectrum) collected within each pixel as a digital number (DN). These DNs represent the surface reflectance from the Earth plus contributions from atmospheric gas absorption, atmospheric scattering, variations in illumination from topographical features, instrument response curves and other artifacts. This is to say, DNs are not surface reflectance! DNs are data that require *radiometric processing* in order to obtain images with physically meaningful quantities like radiance or reflectance. On top of that, depending on the sensor architecture, the pixel values may need to be stitched together to form a full image, and image distortions such as those from sensor tilt or topographical features need to be corrected for. The images also need to be referenced to a real spot on the ground. These corrections are called *geometric corrections*. Hyperspectral imagery, just like any other remote sensing imagery, can be processed or corrected to different levels. These levels are standardized (to some extent) and are used to succinctly describe the stage at which the data is prepared, from raw data from the sensor to mapped variables, and beyond. ## Processing Levels Earth imagery is processed across five main levels, ranging from Level 0 to Level 4. While the precise definitions may differ depending on who you talk to, exactly whose data you're using, or the type of data, they generally follow the high-level framework given below. | Syntax | Description | Wyvern Data | | ----------- | ----------- | :-----------: | | Level 0 (L0) | Reconstructed, unprocessed instrument data at full resolution with all artifacts removed. It should include the image data itself, image metadata, and additional relevant ancillary data with the data . | | | Level 1A (L1A) | Reconstructed, unprocessed instrument data at full resolution, time-referenced, and annotated with ancillary information (camera settings, spacecraft altitude information, platform orbital ephemeris data, etc.). | | | Level 1B (L1B) | Image data which has been radiometrically processed to sensor units (at-sensor radiance, or top-of-atmosphere (TOA) radiance). Includes geolocation reference information and calibration parameters required to reach this level of processing. | :heavy_check_mark: | | Level 2A (L2A) | Image data which has been atmospherically compensated/corrected to bottom-of-atmosphere (BOA) surface reflectance. | :heavy_check_mark: | | Level 2B (L2B) | L2B data are L2A data that have been processed to instrument units. | | | Level 3 (L3) | Data that has been mapped on a uniform space-time grid, and thus have been averaged over space and/or time. | | | Level 4 (L4) | Model outputs or results from analyses of lower level data, potentially from several data sources. | | ## Wyvern's L1B Data Wyvern's hyperspectral imagery data products that are delivered in L1B processing level are both geometrically corrected via georeferencing to map-projected north-up Geographic WGS84 ([EPSG:4326](https://epsg.io/4326)) coordinate system and radiometrically corrected to at-sensor radiance pixel units based on spacecraft location plus pointing along with solar conditions at time of data acquisition. Wyvern's L1B processing level imagery product is delivered as a raster dataset with 32-bit floating point (float32) data type where the pixels represent at-sensor top-of-atmosphere (TOA) radiance that is already in spectral radiance units (W m-2 sr-1 Β΅m-1) with no need to apply any scaling factor. Wyvern's initial L1B imagery product is being delivered with adequate metadata to enable industry standard radiometric & geometric data processing along with advanced analytical techniques such as deep learning, machine learning, spectral indices, classification, anomaly detection, material identification, sub-pixel mixture analysis, and spectral target detection. Since our L1B imagery data product is delivered in Geographic WGS84 (EPSG:4326) coordinate system the spatial X & Y cell size of any given image pixel is defined in angular longitude & latitude degree units. In order to preserve the inherent ground sample distance (GSD) of the data, the pixel cell size for any given imagery data product is set as the longitude & latitude angular degree equivalents of X = 5.0 m x Y = 5.0 m spatial size for the latitude location of the center of the image scene. Consequently, the pixel cell size in longitude & latitude degree units will be non-square and vary based on the latitude location of any given imagery data product (e.g., at 40Β°N latitude X = 5.0 m x Y = 5.0 m in cartesian map distance translates to Lon=0.0000586Β° x Lat=0.0000450Β° in geographic angular degrees). Wyvern's L1B data is delivered in cloud-optimized [GeoTIFF](hyperspectral_image_formats?id=geotiff-tiff) (COG) alongside [STAC](stac_documentation_and_tools) .json metadata. :::info More information about our L1B Data Product can be found in the our [Product Guide](https://guide.wyvern.space/data-product-guide/product-specifications/) ::: ## Using L1B Data Depending on your application, L1B data may still require significant processing to make it suitable for your analysis. Over and above functional requirements, like reprojecting the data to another coordinate reference system (CRS) if it is not compatible with your existing study, or extracting your particular area of interest (AOI), keep in mind that L1B data represents TOA radiance. Absorption of light by ozone, water, carbon dioxide, and other molecules in the atmosphere prevents specific wavelengths from reaching the sensor. Accordingly, the signal-to-noise ratio (SNR) of TOA imagery in wavelengths that overlap with *atmospheric absorption bands* will be lower. This effect is called the [atmospheric window](https://gisgeography.com/atmospheric-window/). L1B compared to L2A calculated via ENVI`s QUAC (QUick Atmospheric Correction) tool. This atmospheric absorption needs to be corrected in order to retrieve surface reflectance values (BOA reflectance). A tool like Py6S or specialized software like ENVI can help you perform atmospheric corrections. We [go into more detail in the atmospheric corrections page](atmospheric_correction). ## 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/index_calculation] # Index Calculations with Python :::warning This tutorial is out of date and uses old hyperspectral drone data from before Wyvern had operating satellites. It still has useful information, but will be replaced soon. ::: 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, however with Python we have to do more of the hard work ourselves! Index calculations, in their simplest form, is just math applied to arrays of numbers (raster images!). There are a number of tools available in Python for working with arrays, and we'll cover the most popular (Numpy) here. :::info Try running this code in a Python notebook to easily troubleshoot & iterate on results! ::: ## Numpy Numpy is fantastic at array calculations, however Numpy is incapable of reading our raster files into arrays. For this, we'll use Rasterio (which relies on GDAL), however there are a number of Python packages that don't have a GDAL dependency that can also be used. :::info We'll be using part of our SFB 15/16 dataset for this analysis. ::: First, we'll read our array into Numpy using Rasterio. ```python import rasterio FILE_PATH = "/vsis3/wyvern-drone-data-processed/processed/SFB1516/WYVERN_SFB1516_20220726_v0_1_0_5mGSD.tiff" file = rasterio.open(FILE_PATH) image_array = file.read() # Leaving band number empty in read operation will read entire raster to array # We should now have an array with 3 dimensions: Band, X, Y (or some combination thereof) print(image_array.shape) ``` Now that we've read our raster into a three dimensional Numpy array, we can easily perform calculations on the array just like we would if calculating a single number. Here we've also defined the array index of each band as a separate variable to make things clearer. We're going to calculate RENDVI, a relatively simple index to calculate. ```python # CHANGE THESE TO THE INDICES OF 750nm AND 705nm BAND IN YOUR FILE band750 = 23 band705 = 19 # RENDVI calculation rendvi_array = (image_array[band750] - image_array[band705]) / (image_array[band750] + image_array[band705]) ``` It's that easy! If we plot the resulting array next to RGB bands, we would see the following image: ![RGB image of a field, next to a green and yellow image of the same field. Green spots indicate healthy vegetation, and yellow areas of less healthy vegetation or dirt](images/rendvi_code_output.png 'Calculated RENDVI vs. RGB :size=800') There is a a more detailed example of index calculations and plotting [available on the Plant Health Indices page.](/docs/documentation/datasets_and_use_cases/plant_health_indices?id=do-it-yourself) ## Numerical Expressions Numpy is fantastic for simple array calculations, however when we look at more complex equations it can get very hard to understand, very quickly. We also would run into issues if we wanted to calculate multiple indices, and automate the process of finding band centers. [Luckily, there is an additional package (NumExpr)](https://github.com/pydata/numexpr) that allows us to define our equations using a string, and provide our band indexes as a dictionary. We can use the same read operation as in the above example: ```python import rasterio FILE_PATH = "/vsis3/wyvern-drone-data-processed/processed/SFB1516/WYVERN_SFB1516_20220726_v0_1_0_5mGSD.tiff" file = rasterio.open(FILE_PATH) image_array = file.read() # Leaving band number empty in read operation will read entire raster to array # We should now have an array with 3 dimensions: Band, X, Y (or some combination thereof) print(image_array.shape) ``` We can then define our equations and band mapping dictionaries: ```python RENDVI_EQUATION = "(b750 - b705) / (b750 + b705)" VREI2_EQUATION = "(b734 - b747) / (b715 + b726)" MTVI2_EQUATION = "(1.5 * (1.2 * (b800 - b550) - 2.5 * (b670 - b550))) / sqrt((2 * b800 + 1) ** 2 - (6 * b800 - 5 * sqrt(b670)) - 0.5)" band_map_dictionary = { "b750": 23, "b705": 19, "b734": 22, "b747": 23, "b715": 20, "b726": 21, "b800": 26, "b550": 8, "b670": 16, } ``` We can then use NumExpr to evaluate these numerical expressions quickly and easily! ```python import numexpr as ne rendvi_array = ne.evaluate(RENDVI_EQUATION, local_dict=band_map_dictionary) vrei2_array = ne.evaluate(VREI2_EQUATION, local_dict=band_map_dictionary) mtvi_array = ne.evaluate(MTVI2_EQUATION, local_dict=band_map_dictionary) ``` Plotting these using indices using matplotlib, we can see that the RENDVI is consistent with what was observed in the above Numpy example. ![Three plots, RENDVI, VREI2, and MTVI2, respectively. RENDVI is mostly green, while VREI2 and MTVI2 show more green around areas of healthy plant growth](images/numexpr.png 'RENVI, VREI2, and MTVI2 plotted :size=1000') ## 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 userful `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/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`. :::info Try running this code in a Python notebook to easily troubleshoot & iterate on results! ::: 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](https://knowledge.wyvern.space/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](https://knowledge.wyvern.space/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 our [hyperspectral library](../../../hyperspectral_library) 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 We leverage the [SpatioTemporal Asset Catalog (STAC)](https://stacspec.org/en/about/) format to make our geospatial data more accessible, interoperable, and easier to discover. ## What is STAC? The SpatioTemporal Asset Catalog (STAC) is an open source standard for organizing and describing geospatial data, ensuring datasets are easily searchable, shareable, and usable. Using a standard such as STAC ensures that you can use multiple datasets within the same analysis easily, whether you’re working with satellite imagery, aerial photos, or other Earth observation data. STAC therefore works to reduce the friction in integrating new geospatial information into your workflows. A *SpatioTemporal Asset* is any file that represents information about the Earth captured in a certain place and at a particular time. For example, that could be any type of imagery (hyperspectral, multispectral, etc.) from satellites, drones, or airplanes, point clouds from LiDAR instruments, or full-motion video. STAC is a network of JavaScript Object Notation (JSON) files that reference other JSON files, with each JSON file adhering to a specific core specification depending on which STAC component it is describing. This core JSON format can also be customized to fit differing needs, making the STAC specification highly flexible and adaptable`[1]`. The JSON files making up the STAC catalog serve like a map to find the STAC Assets (in our case, the hyperspectral images and their metadata). ## Key components of STAC ### Item An Item is a single geospatial asset. It contains metadata such as the time and place the asset was captured and links to the data itself (e.g., GeoTIFF files and thumbnails). A STAC Item represents assets as [GeoJSON](https://en.wikipedia.org/wiki/GeoJSON) with additional metadata so that it can be easily searched. Because an item is formatted as a GeoJSON, it can be easily read by any modern GIS or geospatial library. > **Example Use**: A single satellite image taken on a specific date over a region: "Image of Toronto captured on 2023-03-15." The STAC Item JSON specification employs standard GeoJSON fields along with additional informational fields to provide a more comprehensive description of the asset(s). Figure 1 below shows the fields of the STAC Item JSON. Figure 1: STAC Item Specification. The fields in the gray text (stac_extention and collection) represent optional fields `[2]` ### Catalog The Catalog is the backbone of STAC, organizing geospatial assets into a hierarchy for browsing and discovery. Think of it as the table of contents or a map for all available datasets. > **Example Use**: A Catalog could represent all datasets collected by a satellite mission or stored by a specific provider: "Dragonette-001 Data" or "Wyvern Hyperspectral Data" A catalog.json file contains links to some combination of other STAC Catalogs, Collections, and/or Items. Figure 2 describes the required fields in a STAC Catalog JSON. Figure 2: STAC Catalog Specification. The fields in white text (title and stac_extentions) represent optional fields. `[2]` ### Collection Collections group related Items that share common metadata. For instance, all images from a specific period, sensor, or data acquisition campaign may be part of a Collection. > **Example Use**: Grouping temporal observations from a specific satellite: "Imagery from Dragonette-001 on 2023-03-15" or "Imagery from Dragonette-001 from 2024" A STAC Collection includes all the functionality of a [catalog](#catalog) (organizing and linking items or sub-catalogs), however it adds required metadata fields like extent, license, and providers to better describe the items in the collection. While a STAC Catalog is generic, a STAC Collection is designed to be more informative and specific to its content as it focuses on a coherent group of items with shared characteristics. Figure 3 shows how the STAC Collection specification adds additional (optional) fields which can better describe the group of items. Figure 3: STAC Collection Specification. In addition to the optional fields from the STAC Catalog (title, stac_extensions), the fields providers, keywords, assets, and summaries represent optional fields for Collections. `[2]` ## Why Wyvern uses STAC Wyvern's imagery data products are delivered in a standardized format to make working with the hyperspectral imagery (HSI) data files *as easy as possible*. Wyvern's hyperspectral imagery data products are delivered as Cloud-Optimized GeoTIFF (COG) format raster dataset with embedded image metadata tags along with a separate metadata sidecar file in JSON format that is compliant with the STAC standard. :::info Ready to get started? Check out our [Open Data Portal](https://opendata.wyvern.space/) or [contact us.](mailto:info@wyvern.space) ::: ## Getting Started with STAC You've read this far, so congrats! We'd say you're well on your way to getting started with STAC. Here are some popular tools and libraries you can use to access and interact with STAC-compliant datasets: * [PySTAC](https://pystac.readthedocs.io/en/stable/): A Python library for working with STAC data. It allows you to read, write, and manipulate STAC catalogs, collections, and items. * [STAC Browser](https://github.com/radiantearth/stac-browser): A web-based tool for browsing and searching STAC catalogs. * [sat-utils](https://github.com/sat-utils): A suite of tools and libraries for interacting with satellite data, including STAC-compliant datasets. * [GDAL](https://knowledge.wyvern.space/?_gl=1*1if7u1s*_gcl_au*NTYxMTA5NDEuMTczMzQzMzQ3OQ..#/getting_started/common_tools?id=gdal): While not exclusively STAC-focused, GDAL can interact with geospatial data and can be integrated into workflows using STAC metadata. These tools make it easy to discover, analyze, and integrate STAC data into your workflows. ### 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 Hyperspectral imaging uses a narrow range of wavelengths across the electromagnetic spectrum, to produce images of the Earth's surface. These images contain information about the reflectance of each pixel across the electromagnetic spectrum. How each pixel reflects light can tell us a lot about the materials that are present within a pixel. ## Background We can analyze hyperspectral images through the use of material signatures, also known as spectral signatures or spectral libraries. A material signature is a plot of the reflectance of a particular material as a function of wavelength. These material signatures are typically produced by isolating a specific material, and imaging it under lab conditions. By comparing the reflectance of a pixel in a hyperspectral image to known material signatures, it is possible to identify the materials present within an image. Example of material signatures from the USGS Spectral Library 7 ### Mixed Pixels It is not always possible to accurately identify a material based on the reflectance of a single pixel in a hyperspectral image. This is because a single pixel may contain more than one material, a phenomenon known as a mixed pixel. In these cases, it is necessary to use techniques such as __endmember extraction__ to identify the individual materials present in the mixed pixel. Example of low spatial resolution pixels (black boxes) that contain mixed pixels and unmixed pixels Lower spatial resolutions will result in more mixed pixels as a single pixel covers more of the earth's surface. When spatial resolution is increased, it is more likely that each pixel will contain a single material. Endmember extraction (also known as spectral unmixing) involves identifying a set of pure spectral signatures, known as endmembers, that can be used to represent the materials present in the image. These endmembers are then used to model the mixed pixels in the image, allowing the individual materials within a mixed pixel to be identified and mapped. [PySptools](https://pysptools.sourceforge.io/eea.html) provides endmembers extraction algorithms (EEA) and documentation on how to apply the algorithms. ### Classifying Pixels We can also classify pixels based on known spectral signatures. There are a number of different supervised classification algorithms available, which analyze the similarity of a pixel's reflectance spectra to a known library of material spectral signatures Since most material spectra are collected under laboratory conditions, actual material spectra within hyperspectral imagery may not match exactly. Most classification algorithms for hyperspectral imagery will look at the similarity of an pixel's spectra to known materials, instead of trying to match materials exactly. [PySptools](https://pysptools.sourceforge.io/classification.html) also provides a number of classification algorithms that can be used. --- [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 It's easy to get started in geospatial & hyperspectral processing due to the abundance of free and open source tools available! This is a short list on some of the most common tools available, however the geospatial community is filled with excellent packages and tools for you to explore. ## QGIS Free GUI Open Source [QGIS](https://qgis.org) is a free, open-source geospatial toolkit that is widely used by GIS professionals and researchers for working with spatial data. It allows users to view, edit, and analyze geospatial data, and create custom maps and graphics. QGIS supports a wide range of file formats and is compatible with other open-source and commercial GIS software. It also has a large and active user community that contributes to its development and provides support and resources for users. QGIS is a powerful and versatile tool for working with spatial data, especially considering that it is completely free software. ## ArcGIS Paid GUI Closed Source [ArcGIS](https://www.arcgis.com/home/index.html) is a powerful and widely used GIS software developed by Esri. It allows users to create, manage, and analyze spatial data, and create custom maps and graphics. ArcGIS has a user-friendly interface and a wide range of tools and functions for working with spatial data. It also has strong integration with other Esri products and services, such as ArcGIS Online and ArcGIS Enterprise. ArcGIS is widely used by GIS professionals, businesses, governments, and organizations around the world. It is available in several different versions, including a desktop version and a cloud-based version. ## ENVI Paid GUI Closed Source [ENVI](https://www.nv5geospatialsoftware.com/Products/ENVI) is a geospatial toolkit developed by NV5 Geospatial Solutions for working with hyperspectral imagery. It is widely used by remote sensing professionals, researchers, and organizations for processing, analyzing, and visualizing hyperspectral data. ENVI has a wide range of tools and functions for working with hyperspectral data, including tools for spectral analysis, image processing, and feature extraction. It also has strong integration with other Harris geospatial products, such as IDL and Geomatica. ENVI is available in several different versions, including a standalone version and a version integrated with ArcGIS and other GIS software. ## Python Free Code Open Source [Python](https://www.python.org/) is a popular and versatile programming language that is widely used for geospatial processing and analysis. It has a large and active community of users and developers, and a wide range of tools and libraries specifically designed for working with geospatial data. For example, the Python packages [NumPy](https://numpy.org/) and [Pandas](https://pandas.pydata.org/) are commonly used for working with arrays and tabular data, and the [Python library version of GDAL](https://gdal.org/en/stable/api/python/python_bindings.html) is widely used for working with raster and vector data. Python also has strong integration with other GIS software, such as QGIS and ArcGIS, which allows users to automate and extend the functionality of these tools. Overall, Python is a powerful and versatile tool for working with geospatial data. ## R Free Code Open Source [R](https://www.r-project.org/) is a popular programming language and environment for statistical computing and data analysis, which is widely used by researchers, statisticians, and data scientists. It has an active community of users and developers, and a wide range of tools and packages specifically designed for working with geospatial data. For example, the R packages [sp](https://cran.r-project.org/web/packages/sp/index.html), rgdal, and rgeos are commonly used for working with vector data, and the package [raster](https://cran.r-project.org/web/packages/raster/index.html) is widely used for working with raster data. R also has strong integration with other GIS software, such as QGIS and ArcGIS, which allows users to automate and extend the functionality of these tools. ## GDAL Free Code Open Source [GDAL (Geospatial Data Abstraction Library)](https://gdal.org/) is a widely used open-source library for working with geospatial data. It provides a common interface for working with a variety of different file formats and data sources, and allows users to perform a wide range of geospatial operations, such as projection, resampling, and warping of raster data, and geometry operations on vector data. GDAL is written in C++ and includes bindings for several other programming languages, including Python, R, and C#. It is used by a wide range of GIS software, including QGIS, ArcGIS, and ENVI, and is an essential tool for many geospatial professionals and researchers. --- [Source: https://knowledge.wyvern.space/docs/documentation/getting_started/atmospheric_correction] # Atmospheric Correction In this section, we will talk about two different processing levels, **L1B Top of Atmosphere Radiance**, and **L2A Bottom of Atmosphere Reflectance**. **L1B Top of Atmosphere (ToA) Radiance** is the raw light signal that a satellite sensor measures coming from the Earth, but it still includes the effects of the atmosphere (like scattering and absorption by air and clouds). **L2A Bottom of Atmosphere (BoA) Reflectance**, on the other hand, is a corrected product that shows how much light is actually reflected off the Earth's surface itself, as if the atmosphere wasn't there. This makes L2A easier to use for studying surface features like vegetation, water, and soil, because it better represents what’s really on the ground. Wyvern currently delivers all imagery as L1B Top of Atmosphere, which means that applying atmospheric corrections to the imagery is necessary before any analysis is done. This includes applying indicies like RENDVI. We are currently hard at work building an L2A BoA reflectance product, but in the meantime, there are other tools you can use to atmospherically correct our imagery. The following table lists atmospheric correction tools that are able to easily ingest and process Wyvern data: | Name | Type | Details | | ---- | ---- | ------- | | [ATCOR](https://atcor.com/) | Paid | Trusted atmospheric correction tool. | | [ENVI QUAC](https://www.nv5geospatialsoftware.com/docs/QUAC.html) | Paid | Built-in tool to ENVI. Extremely easy to use. | | [ACOLITE](https://github.com/acolite/acolite) | Free | Open source atmospheric correction tool for coastal applications. | There are also many other tools than can be used with Wyvern data, however manual work would need to be completed to ensure compatibility and correct output. These tools 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), [libRadtran](https://www.libradtran.org/doku.php?id=start), and many more. ## More Resources NASA has [fantastic documentation on processinglevels](https://www.earthdata.nasa.gov/learn/earth-observation-data-basics/data-processing-levels#:~:text=NASA%20data%20products%20are%20processed,data%20at%20full%20instrument%20resolution). --- [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; an [L1B collection](https://wyvern-odp.com/l1b/collection.json) exists but returned zero items as of 2026-08 β€” re-check rather than assuming it is still empty. **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); GSD β‰ˆ 5.2 m - 23 bands (Standard VNIR, Dragonette-1, ~503–799 nm) or 31 bands (Extended VNIR, Dragonette-2/3/4, ~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 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.