RemoteSensingToolbox
RemoteSensingToolbox is a pure Julia package built on top of Rasters.jl for reading, visualizing, and processing remotely sensed imagery. Users may refer to the Tutorials section for examples on how to use this package.
Installation
To install this package, first start the Julia REPL and then open the package manager by typing ]
. You can then download RemoteSensingToolbox
directly from the official Julia repository like so:
(@v1.9) pkg> add RemoteSensingToolbox
Once RemoteSensingToolbox
has been installed, you can import it like any other Julia package. Please note that many features require you to also import the Rasters
package.
using RemoteSensingToolbox, Rasters
Features
This package is a work in progress, which means that new features are being added and existing features are subject to change. To contribute, please create an issue on GitHub or open a pull request. A summary of both existing and planned features is provided below:
Feature | Description | Implemented |
---|---|---|
Reading and Writing | Read layers from a scene and the write results to disk | Yes |
Visualization | Visualize images with various band composites | Yes |
Land Cover Indices | Calculate indices such as MNDWI and NDVI | Yes |
QA and SCL Decoding | Decode Quality Assurance and Scene Classification masks | Yes |
Pixel Masking | Mask pixels to remove objects such as clouds or shadows | Yes |
PCA | Perform PCA analysis, transformation, and reconstruction | Yes |
MNF | Minimum Noise Fraction transformation and reconstruction | Yes |
Signature Analysis | Visualize spectral signatures for different land cover types | Yes |
Land Cover Classification | Exposes an MLJ interface for classifying land cover types | No |
Endmember Extraction | Extract spectral endmembers from an image | No |
Spectral Unmixing | Perform spectral unmixing under a given endmember library | No |
Rasters.jl
RemoteSensingToolbox
is intended to be used in conjunction with the wider Julia ecosystem and as such, seeks to avoid duplicating functinalities provided by other packages. As the majority of methods accept and return AbstractRaster
or AbstractRasterStack
objects, users should be able to call methods from Rasters.jl at any point in the processing pipeline. A summary of common functionalities offered by Rasters
is provided below:
Method | Description |
---|---|
mosaic | Join rasters covering different extents into a single array or file. |
crop | Shrink objects to specific dimension sizes or the extent of another object. |
extend | Extend objects to specific dimension sizes or the extent of another object. |
trim | Trims areas of missing values for arrays and across stack layers. |
resample | Resample data to a different size and projection, or snap to another object. |
mask | Mask a raster by a polygon or the non-missing values of another Raster. |
replace_missing | Replace all missing values in a raster and update missingval. |
extract | Extract raster values from points or geometries. |
zonal | Calculate zonal statistics for a raster masked by geometries. |
Satellites
When working with remotely sensed images, the user is often required to adapt their approach according to the type of sensor that produced it. For example, decoding digital numbers into radiance, reflectance, or temperature requires knowledge about the encoding scheme used by the satellite in question. To address these issues, we provide the AbstractSatellite
type, which encodes various sensor-specific parameters at the type level for several popular platforms. This approach allows the end-user to write generic code that automatically adapts to the sensor being used. For example, if we have bound the variable src
to an instance of Landsat8
, then we can load a cloud mask from the included QA file by calling Raster(src, :clouds)
, or calculate the MNDWI index by calling mndwi(src)
.
SatelliteDataSources.AbstractSatellite
— TypeThe super-type of all satellites.
Sub-types must implement the AbstractSatellite
interface.
SatelliteDataSources.Landsat7
— TypeImplements the AbstractSatellite
interface for Landsat 7.
Supported Bands: :B1
, :B2
, :B3
, :B4
, :B5
, :B7
, :thermal
, :panchromatic
Supported Colors: :blue
, :green
, :red
, :nir
, :swir1
, :swir2
Supported Masks: :dilated_clouds
, :clouds
, :cloud_shadow
, :snow
, :water
SatelliteDataSources.Landsat8
— TypeImplements the AbstractSatellite
interface for Landsat 8.
Supported Bands: :B1
, :B2
, :B3
, :B4
, :B5
, :B6
, :B7
, :thermal1
, :thermal2
, :panchromatic
Supported Colors: :blue
, :green
, :red
, :nir
, :swir1
, :swir2
Supported Masks: :dilated_clouds
, :clouds
, :cloud_shadow
, :snow
, :water
SatelliteDataSources.Landsat9
— TypeImplements the AbstractSatellite
interface for Landsat 9.
Supported Layers: :B1
, :B2
, :B3
, :B4
, :B5
, :B6
, :B7
, :thermal1
, :thermal2
, :panchromatic
Supported Colors: :blue
, :green
, :red
, :nir
, :swir1
, :swir2
Supported Masks: :dilated_clouds
, :clouds
, :cloud_shadow
, :snow
, :water
SatelliteDataSources.Sentinel2
— TypeImplements the AbstractSatellite
interface for Sentinel 2.
The user must specify a resolution of 10, 20, or 60 meters.
Supported Bands (10m): :B02
, :B03
, :B04
, :B08
Supported Bands (20m): :B02
, :B03
, :B04
, :B05
, :B06
, :B07
, :B8A
, :B11
, :B12
Supported Bands (60m): :B01
, :B02
, :B03
, :B04
, :B05
, :B06
, :B07
, :B8A
, :B09
, :B11
, :B12
Supported Colors (10m): :blue
, :green
, :red
, :nir
Supported Colors (20m and 60m): :blue
, :green
, :red
, :nir
, :swir1
, :swir2
Supported Masks (20m and 60m): :cloud_shadow
, :clouds_med
, :clouds_high
, :cirrus
, :vegetation
, :soil
, :water
, :snow
SatelliteDataSources.DESIS
— TypeImplements the AbstractSatellite
interface for DESIS.
Supported Bands: :Bands
, :Band_30
, :Band_65
, :Band_100
, :Band_175
Supported Colors: :blue
, :green
, :red
, :nir
Supported Masks: :clouds
, :shadow
, :haze
, :snow
, :land
, :water
SatelliteDataSources.bands
— Functionbands(::Type{AbstractSatellite})
Return the band names in order from shortest to longest wavelength.
SatelliteDataSources.layers
— Functionlayers(::Type{AbstractSatellite})
layers(x::AbstractSatellite)
Return the names of all layers available for the given sensor.
Example
# Get all Available Layers for Landsat 8
landsat_layers = layers(Landsat8)
# Get all Available Layers for a Specific Scene
src = Landsat8("LC08_L2SP_043024_20200802_20200914_02_T1")
available_layers = layers(src)
SatelliteDataSources.wavelengths
— Functionwavelengths(::Type{AbstractSatellite})
Return the central wavelengths for all bands in order from shortest to longest.
SatelliteDataSources.wavelength
— Functionwavelength(::Type{AbstractSatellite}, band::Symbol)
Return the central wavelength for the corresponding band.
SatelliteDataSources.blue_band
— Functionblue_band(::Type{AbstractSatellite})
Return the blue band for the given sensor.
SatelliteDataSources.green_band
— Functiongreen_band(::Type{AbstractSatellite})
Return the green band for the given sensor.
SatelliteDataSources.red_band
— Functionred_band(::Type{AbstractSatellite})
Return the red band for the given sensor.
SatelliteDataSources.nir_band
— Functionnir_band(::Type{AbstractSatellite})
Return the nir band for the given sensor.
SatelliteDataSources.swir1_band
— Functionswir1_band(::Type{AbstractSatellite})
Return the swir1 band for the given sensor.
SatelliteDataSources.swir2_band
— Functionswir2_band(::Type{AbstractSatellite})
Return the swir2 band for the given sensor.
SatelliteDataSources.dn_scale
— Functiondn_scale(::Type{AbstractSatellite}, layer::Symbol)
Return the scale factor applied to digital numbers.
SatelliteDataSources.dn_offset
— Functiondn_offset(::Type{AbstractSatellite}, layer::Symbol)
Return the offset factor applied to digital numbers.
SatelliteDataSources.decode
— Functiondecode(s::Type{AbstractSatellite}, raster::Rasters.AbstractRaster)
decode(s::Type{AbstractSatellite}, raster::Rasters.AbstractRasterStack)
Decode the Digital Number (DN) values in the given raster(s). Typically, the decoded values will be in either reflectance (visual bands) or Kelvin (thermal bands).
Parameters
s
: TheAbstractSatellite
that produced the raster(s) in question.raster
: Either aRasters.Raster
orRasters.RasterStack
to be decoded.
SatelliteDataSources.encode
— Functionencode(s::Type{AbstractSatellite}, raster::Rasters.AbstractRaster; encoding_type=UInt16, missingval=0x0000)
encode(s::Type{AbstractSatellite}, raster::Rasters.AbstractRasterStack; kwargs...)
Encode the provided raster(s) to Digital Numbers (DN).
Parameters
s
: TheAbstractSatellite
that produced the raster(s) in question.raster
: Either aRasters.Raster
orRasters.RasterStack
to be encoded.encoding_type
: The Julia type to use for storing the encoded values (default =UInt16
).missingval
: The value to denote missing values (default =0x0000
).
SatelliteDataSources.metadata
— Functionmetadata(x::AbstractSatellite)
Parses the metadata fields for the given satellite scene.
Metadata varies between products, but typically includes the acquisition date and processing level.
Rasters.Raster
— TypeRaster(x::AbstractSatellite, layer::Symbol; kwargs...)
Read a single layer into a Rasters.Raster
.
Parameters
x
: AnAbstractSatellite
from which to read a layer.layer
: The layer to be read. Seelayers(s)
for a list of available layers for sensors
.kwargs
: Refer to theRasters.Raster
documentation for a summary of supported keywords.
Rasters.RasterStack
— TypeRasterStack(x::AbstractSatellite, layers=bands(T); kwargs...)
Read multiple layers into a Rasters.RasterStack
.
Parameters
x
: AnAbstractSatellite
from which to read a layer.layer
: The layer to be read. Seelayers(s)
for a list of available layers for sensors
.kwargs
: Refer to theRasters.RasterStack
documentation for a summary of supported keywords.
Visualization
Remotely sensed imagery is typically encoded as either UInt16
or Int16
values. However, many products only actually use the first 12 bits for storing information. The result is that naive visualization methods will produce a near-black image, since the maximum possible brightness will be located in the lower range of values provided by the 16 bit encoding. To address this, we need to perform a linear stretch before visualizing an image. Moreover, many satellites have more than three bands, which motivates the use of band combinations to emphasize certain features and land cover types.
RemoteSensingToolbox.visualize
— Functionvisualize(g::AbstractRaster; lower=0.02, upper=0.98)
visualize(r::AbstractRaster, g::AbstractRaster, b::AbstractRaster; lower=0.02, upper=0.98)
Visualize an RGB or grayscale satellite image.
Returns an Array
of either RGB{N0f8}
or Gray{N0f8}
elements.
Keywords
lower
: The lower percentile to use for adjusting the image histogram.upper
: The upper percentile to use for adjusting the image histogram.
Example
using RemoteSensingToolbox, Rasters, FileIO
# Prepare a Landsat 8 Image
src = Landsat8("LC08_L2SP_043024_20200802_20200914_02_T1")
# Display True Color Image
r, g, b = RasterStack(src, [:red, :green, :blue])
img = visualize(r, g, b; upper=0.90)
FileIO.save("truecolor.jpg", img)
RemoteSensingToolbox.true_color
— Functiontrue_color(src::AbstractSatellite; lower=0.02, upper=0.98)
true_color(s::Type{AbstractSatellite}, raster::AbstractRasterStack; lower=0.02, upper=0.98)
Visualize a satellite image using the true-color band combination, which appears like a natural image.
Accepts either an AbstractSatellite
or a combination of Type{AbstractSatellite}
and AbstractRasterStack
.
Returns an Array
of either RGB{N0f8}
or Gray{N0f8}
elements.
Keywords
lower
: The lower percentile to use for adjusting the image histogram.upper
: The upper percentile to use for adjusting the image histogram.
Example
using RemoteSensingToolbox, Rasters
# Prepare a Landsat 8 Image
src = Landsat8("LC08_L2SP_043024_20200802_20200914_02_T1")
# Read Bands Directly from Disk
true_color(src; upper=0.90)
# Read Bands from a RasterStack
stack = RasterStack(src; lazy=true)
true_color(Landsat8, stack; upper=0.90)
RemoteSensingToolbox.color_infrared
— Functioncolor_infrared(src::AbstractSatellite; lower=0.02, upper=0.98)
color_infrared(s::Type{AbstractSatellite}, raster::AbstractRasterStack; lower=0.02, upper=0.98)
Visualize a satellite image using the color-infrared band combination, which highlight vegetation in red, water in blue, and urban areas in grey.
Accepts either an AbstractSatellite
or a combination of Type{AbstractSatellite}
and AbstractRasterStack
.
Returns an Array
of either RGB{N0f8}
or Gray{N0f8}
elements.
Keywords
lower
: The lower percentile to use for adjusting the image histogram.upper
: The upper percentile to use for adjusting the image histogram.
RemoteSensingToolbox.swir
— Functionswir(src::AbstractSatellite; lower=0.02, upper=0.98)
swir(s::Type{AbstractSatellite}, raster::AbstractRasterStack; lower=0.02, upper=0.98)
Visualize a satellite image using the SWIR band combination, which emphasizes dense vegetation as dark green.
Accepts either an AbstractSatellite
or a combination of Type{AbstractSatellite}
and AbstractRasterStack
.
Returns an Array
of either RGB{N0f8}
or Gray{N0f8}
elements.
Keywords
lower
: The lower percentile to use for adjusting the image histogram.upper
: The upper percentile to use for adjusting the image histogram.
RemoteSensingToolbox.agriculture
— Functionagriculture(src::AbstractSatellite; lower=0.02, upper=0.98)
agriculture(s::Type{AbstractSatellite}, raster::AbstractRasterStack; lower=0.02, upper=0.98)
Visualize a satellite image with the agricultural band combination, which is used for crop monitoring and emphasizes healthy vegetation.
Accepts either an AbstractSatellite
or a combination of Type{AbstractSatellite}
and AbstractRasterStack
.
Returns an Array
of either RGB{N0f8}
or Gray{N0f8}
elements.
Keywords
lower
: The lower percentile to use for adjusting the image histogram.upper
: The upper percentile to use for adjusting the image histogram.
RemoteSensingToolbox.geology
— Functiongeology(src::AbstractSatellite; lower=0.02, upper=0.98)
geology(s::Type{AbstractSatellite}, raster::AbstractRasterStack; lower=0.02, upper=0.98)
Visualize a satellite image with the geology band combination, which emphasizes geological formations, lithology features, and faults.
Accepts either an AbstractSatellite
or a combination of Type{AbstractSatellite}
and AbstractRasterStack
.
Returns an Array
of either RGB{N0f8}
or Gray{N0f8}
elements.
Keywords
lower
: The lower percentile to use for adjusting the image histogram.upper
: The upper percentile to use for adjusting the image histogram.
Land Cover Indices
Land cover indices are used to highlight different types of land cover. For example, the Modified Normalized Difference Water Index (MNDWI) is used to highlight water while diminishing built-up areas. Each index is expressed as a function of two or more bands. RemoteSensingToolbox
can automatically select the appropriate bands for a given index by specifying the AbstractSatellite
type, while also providing the option to manually specify bands as the user desires.
RemoteSensingToolbox.mndwi
— Methodmndwi(src::AbstractSatellite)
mndwi(::Type{AbstractSatellite}, stack::AbstractRasterStack)
mndwi(green::AbstractRaster, swir::AbstractRaster)
Compute the Modified Normalised Difference Water Index (Xu 2006).
MNDWI = (green - swir) / (green + swir)
RemoteSensingToolbox.nbri
— Methodnbri(src::AbstractSatellite)
nbri(::Type{AbstractSatellite}, stack::AbstractRasterStack)
nbri(nir::AbstractRaster, swir2::AbstractRaster)
Compute the Normalized Burn Ratio Index.
NBRI is used to emphasize burned areas.
NBRI = (nir - swir2) / (nir + swir2)
RemoteSensingToolbox.ndbi
— Methodndbi(src::AbstractSatellite)
ndbi(::Type{AbstractSatellite}, stack::AbstractRasterStack)
ndbi(swir1::AbstractRaster, nir::AbstractRaster)
Compute the The Normalized Difference Built-up Index
NDBI is used to emphasize urban and built-up areas.
NDBI = (swir1 - nir) / (swir1 + nir)
RemoteSensingToolbox.ndmi
— Methodndmi(src::AbstractSatellite)
ndmi(::Type{AbstractSatellite}, stack::AbstractRasterStack)
ndmi(nir::AbstractRaster, swir1::AbstractRaster)
Compute the Normalized Difference Moisture Index.
NDMI is sensitive to the moisture levels in vegetation. It is used to monitor droughts and fuel levels in fire-prone areas.
NDMI = (nir - swir1) / (nir + swir1)
RemoteSensingToolbox.ndvi
— Methodndvi(src::AbstractSatellite)
ndvi(::Type{AbstractSatellite}, stack::AbstractRasterStack)
ndvi(nir::AbstractRaster, red::AbstractRaster)
Compute the Normalized Difference Vegetation Index.
NDVI = (nir - red) / (nir + red)
RemoteSensingToolbox.ndwi
— Methodndwi(src::AbstractSatellite)
ndwi(::Type{AbstractSatellite}, stack::AbstractRasterStack)
ndwi(green::AbstractRaster, nir::AbstractRaster)
Compute the Normalized Difference Water Index (McFeeters 1996).
NDWI = (green - nir) / (green + nir)
RemoteSensingToolbox.savi
— Methodsavi(src::AbstractSatellite; L=0.33)
savi(::Type{AbstractSatellite}, stack::AbstractRasterStack; L=0.33)
savi(nir::AbstractRaster, red::AbstractRaster; L=0.33, scale=1.0, offset=0.0)
Compute the Soil Adjusted Vegetation Index (Huete 1988).
SAVI is a vegetation index which attempts to minimize soil brightness influences by introducing a soil-brightness correction factor (L).
SAVI = ((nir - red) / (nir + red + L)) * (1 + L)
Keywords
L
: The ammount of vegetative cover, where 1.0 means no vegetation and 0.0 means high vegetation.scale
: The scaling factor to convert digital numbers to reflectance.offset
: The offset to convert digital numbers to reflectance.
Spectral Analysis
Spectral analysis involves studying the relationships between different materials and their corresponding spectral signatures. Due to the interactions between light and matter, each signature is unique to the material that emitted it. We can exploit this fact to assign a label to each pixel, or even estimate the abundances of different materials at a sub-pixel level. This package provides several methods for both extracting and visualizing these signatures.
RemoteSensingToolbox.extract_signatures
— Functionfunction extract_signatures([agg], raster, shp, label::Symbol)
Extract signatures from a Raster
or RasterStack
within regions specified by a provided shapefile.
Parameters
agg
: A function to aggregate signatures belonging to the same class (ex.mean
,median
,maximum
).raster
: AnAbstractRaster
orAbstractRasterStack
from which to extract the spectral signatures.shp
: A table with a:geometry
column ofGeoInterface.jl
geometries and land cover labels.label
: The column inshp
corresponding to the land cover type.
Returns
A Tables.columntable
containing all labelled signatures or their aggregation if agg
is provided.
Example
julia> src = Landsat8("LC08_L2SP_043024_20200802_20200914_02_T1/");
julia> sr = decode(Landsat8, RasterStack(src, [:B1, :B2, :B3, :B4, :B5]));
julia> shp = GeoDataFrames.read("data/landcover/landcover.shp");
julia> extract_signatures(sr, shp, :C_name) |> DataFrame
1925×6 DataFrame
Row │ label B1 B2 B3 B4 B5
│ String Float32 Float32 Float32 Float32 Float32
──────┼───────────────────────────────────────────────────────────────
1 │ Hail Scar 0.058005 0.10613 0.19028 0.25177 0.51577
2 │ Hail Scar 0.057895 0.10888 0.19358 0.257215 0.520198
3 │ Hail Scar 0.06026 0.111795 0.197265 0.263045 0.523855
4 │ Hail Scar 0.0595175 0.109375 0.195257 0.258315 0.522287
5 │ Hail Scar 0.059215 0.108468 0.191655 0.254327 0.51401
⋮ │ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮
1921 │ Built Up 0.0633125 0.0914725 0.145565 0.165227 0.255153
1922 │ Built Up 0.0764025 0.100878 0.154392 0.176063 0.24363
1923 │ Built Up 0.09788 0.126342 0.180738 0.19941 0.25958
1924 │ Built Up 0.0973025 0.133163 0.1894 0.198475 0.278775
1925 │ Built Up 0.06345 0.089685 0.145042 0.158243 0.272505
1915 rows omitted
julia> extract_signatures(mean, sr, shp, :C_name) |> DataFrame
7×6 DataFrame
Row │ label B1 B2 B3 B4 B5
│ String Float32 Float32 Float32 Float32 Float32
─────┼─────────────────────────────────────────────────────────────────────
1 │ Hail Scar 0.0618124 0.107894 0.187857 0.24683 0.507976
2 │ Bare Earth 0.0484324 0.053983 0.0775185 0.0886341 0.231951
3 │ Road 0.0400183 0.0534427 0.0938634 0.0958244 0.245025
4 │ Lake 0.00137886 0.00423271 0.0135606 0.00652965 0.031825
5 │ Trees 0.0183338 0.0203976 0.0401544 0.024148 0.318622
6 │ Vegetation 0.00441971 0.0157759 0.0787841 0.0463343 0.495117
7 │ Built Up 0.0845031 0.113944 0.174032 0.197151 0.283717
RemoteSensingToolbox.plot_signatures
— Functionplot_signatures(bandset::Type{<:AbstractSatellite}, sigs; kwargs...)
plot_signatures(bandset::Vector{<:Pair}, sigs; colors=Makie.wong_colors(), label=:label)
Plot the spectral signatures for one or more land cover types.
Parameters
bandset
: AnAbstractSatellite
or a vector of sortedband => wavelength
pairs.sigs
: A table whose rows consist of spectral signatures and their associated labels.
Keywords
label
: The column insigs
containing the signature labels (default =:label
).colors
: The color scheme used by the plot (default =Makie.wong_colors()
).
Example
julia> src = Landsat8("LC08_L2SP_043024_20200802_20200914_02_T1");
julia> surface_reflectance = decode(Landsat8, RasterStack(src));
julia> shp = GeoDataFrames.read("data/landcover/landcover.shp");
julia> sigs = extract_signatures(mean, surface_reflectance, shp, :MC_name);
julia> plot_signatures(Landsat8, sigs, colors=[:brown, :orange, :blue, :green])
RemoteSensingToolbox.plot_signatures!
— Functionplot_signatures!(ax, bandset::Type{<:AbstractSatellite}, sigs; kwargs...)
plot_signatures!(ax, bandset::Vector{<:Pair}, sigs; colors=Makie.wong_colors(), label=:label)
Plot the spectral signatures for one or more land cover types onto an existing Makie.Axis
.
Parameters
ax
: TheMakie.Axis
onto which to plot the signatures.bandset
: AnAbstractSatellite
or a vector of sortedband => wavelength
pairs.sigs
: A table whose rows consist of spectral signatures and their associated labels.
Keywords
label
: The column insigs
containing the signature labels (default =:label
).colors
: The color scheme used by the plot (default =Makie.wong_colors()
).
Example
using RemoteSensingToolbox, Rasters, GeoDataFrames, Statistics, CairoMakie
# Read Images And Convert DNs To Reflectance
landsat = Landsat8("data/LC08_L2SP_043024_20200802_20200914_02_T1/")
sentinel = Sentinel2{20}("data/T11UPT_20200804T183919/")
landsat_reflectance = decode(Landsat8, RasterStack(landsat))
sentinel_reflectance = decode(Sentinel2, RasterStack(sentinel))
# Read Landcover Labels From Shapefile
shp = GeoDataFrames.read("data/landcover/landcover.shp")
# Extract Average Spectral Signature For Each Landcover Type
landsat_sigs = extract_signatures(mean, landsat_reflectance, shp, :MC_name)
sentinel_sigs = extract_signatures(mean, sentinel_reflectance, shp, :MC_name)
# Create Figure and Axis
fig = Figure(resolution=(800,550));
ax1 = Axis(fig[1,1], ylabel="Reflectance", title="Landsat Signatures")
ax2 = Axis(fig[2,1], xlabel="Wavelength (nm)", ylabel="Reflectance", title="Sentinel Signatures")
# Plot Signatures
plot_signatures!(ax1, Landsat8, landsat_sigs, colors=[:brown, :orange, :blue, :green])
plot_signatures!(ax2, Sentinel2{20}, sentinel_sigs, colors=[:brown, :orange, :blue, :green])
# Add Legend
Legend(fig[:,2], ax1)
Principal Component Analysis
Principal Component Analysis (PCA) is typically used to reduce the dimensionality of data. In the case of remote sensing, we are interested in reducing the number of bands we need to store while retaining as much information as possible. PCA rotates the bands into a new coordinate space where each band, called a principal component, is orthogonal to and uncorrelated with every other component. By convention, we order the bands in the transformed image in terms of their explained variance, such that the nth component accounts for more variance than any component after it.
RemoteSensingToolbox.PCA
— TypeRemotely sensed imagery typically consists of anywhere from four to several hundred spectral bands. These bands are often highly correlated due to occupying similar spectral regions. Principal Component Analysis (PCA) is used in remote sensing to:
Create a smaller dataset from multiple bands, while retaining as much of the original spectral information as possible. The new image will consist of several uncorrelated PC bands.
Reveal complex relationships among spectral features.
Distinguish between characteristics that are prevalent in most bands and those that are specific to only a few.
Example
julia> desis = DESIS("DESIS-HSI-L2A-DT0884573241_001-20200601T234520-V0210/");
julia> desis_bands = Raster(desis, :Bands)
julia> pca = fit_pca(desis_bands, stats_fraction=0.1)
PCA(dimensions=235)
Projection Matrix:
235×235 Matrix{Float32}:
0.0001 0.0032 0.0016 -0.0094 0.0147 -0.0151 -0.0049 0.0163 … 0.0038 -0.0012 0.0008 -0.0032 -0.0007 0.0053
0.0005 0.0099 0.0042 -0.0244 0.0517 -0.0335 -0.0185 0.0441 -0.0023 -0.0016 -0.0008 0.001 0.0003 -0.0004
0.0003 0.015 0.0053 -0.037 0.133 -0.0443 -0.0271 0.1381 -0.0006 0.0004 0.0002 -0.0006 0.0003 0.0002
0.0003 0.019 0.0071 -0.0385 0.1369 -0.0393 -0.0148 0.0949 -0.0008 0.0006 -0.0001 -0.0006 0.0001 -0.0006
0.0003 0.0232 0.0073 -0.0419 0.1469 -0.037 0.0041 0.0839 -0.0013 -0.0025 0.0005 0.0007 0.0007 -0.0002
0.0001 0.0267 0.0077 -0.0461 0.1713 -0.0325 0.0246 0.0861 … 0.0013 -0.0007 -0.0007 0.0024 -0.0012 -0.0028
0.0001 0.0295 0.0083 -0.0476 0.1695 -0.0348 0.0319 0.0827 -0.0015 -0.0016 -0.0029 0.0004 0.0019 -0.0009
-0.0001 0.0318 0.0086 -0.0482 0.17 -0.0352 0.0414 0.0784 0.0005 -0.002 0.0003 0.0021 -0.0003 -0.0022
⋮ ⋮ ⋱ ⋮
-0.0663 0.0371 -0.1728 0.0196 -0.0508 -0.1394 -0.0054 -0.0226 -0.0003 0.0001 0.0007 -0.0009 -0.0002 0.0004
-0.0658 0.0365 -0.1679 0.0204 -0.0717 -0.1474 -0.0087 -0.0193 -0.0006 0.0002 0.0004 -0.0 0.0004 -0.0001
-0.0655 0.0352 -0.163 0.0193 -0.0767 -0.1511 -0.012 -0.0232 -0.0005 0.0002 -0.0004 -0.0001 -0.0002 0.0005
-0.066 0.035 -0.1618 0.0193 -0.0859 -0.1503 -0.0055 -0.0177 … 0.0003 0.0005 0.001 -0.0 0.0003 -0.0002
-0.067 0.035 -0.1619 0.019 -0.0745 -0.1466 0.0245 -0.0228 -0.0001 0.0002 -0.0002 -0.0001 -0.0002 0.0005
-0.0679 0.0343 -0.1601 0.0176 -0.0721 -0.139 0.0328 -0.0286 0.0003 0.0003 0.0005 -0.0 0.0003 0.0002
-0.0682 0.0337 -0.1588 0.0151 -0.0458 -0.1242 0.0549 -0.0315 0.0002 -0.0002 -0.0003 0.0002 -0.0003 -0.0005
-0.0711 0.0343 -0.1601 0.012 -0.0612 -0.1804 0.2151 -0.0971 0.0004 -0.0001 0.0 0.0003 0.0 -0.0003
Importance of Components:
Cumulative Variance: 0.8782 0.9493 0.9809 0.9869 0.9889 0.9906 0.9915 0.9922 0.9928 ... 1.0
Explained Variance: 0.8782 0.0711 0.0316 0.006 0.0021 0.0017 0.0009 0.0007 0.0006 ... 0.0
julia> transformed = forward_pca(pca, desis_bands, 12);
julia> size(transformed)
(1131, 1120, 12)
julia> recovered = inverse_pca(pca, transformed);
julia> size(recovered)
(1131, 1120, 235)
RemoteSensingToolbox.cumulative_variance
— Methodcumulative_variance(x::PCA)
Return the cumulative variance associated with each principal component of the fitted PCA transform.
RemoteSensingToolbox.explained_variance
— Methodexplained_variance(x::PCA)
Return the explained variance associated with each principal component of the fitted PCA transform.
RemoteSensingToolbox.fit_pca
— Methodfit_pca(raster::Union{<:AbstractRasterStack, <:AbstractRaster}; method=:cov, stats_fraction=1.0)
fit_pca(signatures::Matrix; method=:cov, stats_fraction=1.0)
Fit a Principal Component Analysis (PCA) transformation to the given raster or spectral signatures.
Parameters
raster
: TheAbstractRaster
orAbstractRasterStack
on which to fit the PCA transformation.signatures
: An nxb matrix of spectral signatures where n is the number of signatures and b is the number of bands.
Keywords
method
: Either:cov
or:cor
, depending on whether we want to use the covariance or correlation matrix for computing the PCA rotation.stats_fraction
: The fraction of pixels to use when computing the covariance (or correlation) matrix. Values less than 1.0 will speed up computation at the cost of precision.
RemoteSensingToolbox.forward_pca
— Methodforward_pca(transformation::PCA, raster::Union{<:AbstractRaster, <:AbstractRasterStack}, components::Int)
forward_pca(transformation::PCA, signatures::Matrix, components::Int)
Perform a forward Principal Component Analysis (PCA) rotation while retaining only the specified number of components.
Parameters
transformation
: A previously fitted PCA transformation.raster
: TheAbstractRaster
orAbstractRasterStack
on which to perform the PCA transformation.signatures
: An nxb matrix of signatures where n is the number of signatures and b is the number of bands.components
: The number of bands to retain in the transformed image or signatures.
RemoteSensingToolbox.inverse_pca
— Methodinverse_pca(transformation::PCA, raster::AbstractRaster)
inverse_pca(transformation::PCA, signatures::Matrix)
Perform an inverse Principal Component Analysis (PCA) transformation to recover the original image.
Parameters
transformation
: A previously fitted PCA transformation.raster
: AnAbstractRaster
representing a previously transformed image.signatures
: An nxc matrix of previously transformed signatures where n is the number of signatures and c is the number of components.
RemoteSensingToolbox.projection
— Methodprojection(x::PCA)
Return the projection matrix for the fitted PCA transformation.
Minimum Noise Fraction
The Minimum Noise Fraction (MNF) transformation is used to separate noise from data along the spectral dimension. This method is typically used with hyperspectral imagery, both as a means of dimension reduction and for noise removal. The transformed image will have its bands placed in descending order according to their Signal to Noise Ratio (SNR). The result is that the noise becomes concentrated in the higher bands, which can then be removed by either applying a standard image denoising algorithm or dropping them altogether.
RemoteSensingToolbox.MNF
— TypeThe Minimum Noise Fraction (MNF) transform is a linear transformation used to reduce the spectral dimensionality of image data and segregate noise. MNF consists of two separate principal component rotations. The first rotation uses the principal components of the noise covariance matrix to decorrelate and rescale the noise (a process known as noise whitening), resulting in a transformed image in which the noise has unit variance and no band-to-band correlations. The second rotation is a standard PCA rotation applied to the noise-whitened image.
The bands in the transformed image will be ordered according to their Signal to Noise Ratio (SNR), with the highest SNR being placed first. The result is that the noise becomes concentrated in the higher bands. Thus, the transform can be used to separate noise from data by performing a forward transform, determining which bands contain coherent images, and running an inverse transform after either discarding or denoising the remaining bands. The number of bands to keep in the inverse transform can be determined by a number of methods. The simplest approach is to look at the sorted transformed bands and threshold at the band where no recognizable features can be observed. An alternative method is to threshold at a desired cumulative SNR.
Example
julia> src = DESIS("data/DESIS-HSI-L2A-DT0485529167_001-20220712T223540-V0220")
julia> desis = decode(DESIS, Raster(src, :Bands))
julia> roi = @view desis[X(1019:1040), Y(550:590)];
julia> mnf = fit_mnf(desis, roi)
MNF(dimensions=235)
Projection Matrix:
235×235 Matrix{Float32}:
-2.135 2.3502 0.3612 0.5912 0.5217 -0.0917 0.0043 … 0.0002 -0.0001 -0.0004 -0.0004 -0.0 0.0004
-0.0959 0.0422 -0.0047 -0.2362 -0.3962 -0.2313 -0.1685 0.0001 0.0004 0.0002 0.0003 0.0001 -0.0001
0.0043 0.0058 -0.0032 0.0023 -0.0061 -0.0048 0.0028 -0.0001 0.0001 -0.0001 -0.0 0.0 0.0001
0.0039 0.002 -0.002 0.0012 -0.0032 -0.004 0.006 -0.0006 0.0002 -0.0001 -0.0 -0.0003 -0.0
0.0024 0.0018 -0.003 -0.0008 0.0038 -0.0003 0.0057 0.0009 -0.0007 0.0002 -0.0012 0.0012 -0.0
0.0019 -0.003 0.0038 -0.002 -0.0001 -0.0003 -0.0 … 0.0021 0.0009 -0.0004 0.0014 -0.0011 -0.0006
0.0047 0.0055 0.006 -0.0014 -0.0011 0.0021 0.0053 -0.0035 0.0006 0.0002 -0.0026 0.0014 0.0019
0.0072 0.0042 0.0012 0.0016 0.0011 -0.002 -0.001 0.0 -0.0026 0.0012 0.0034 -0.0006 -0.0009
⋮ ⋮ ⋱ ⋮
-0.0004 -0.0012 0.0007 0.0006 -0.0 0.0002 0.0005 0.0005 0.0002 -0.0004 -0.0002 -0.0001 0.0001
-0.0014 -0.0005 -0.0005 0.0019 -0.0002 -0.0005 0.0017 0.0003 0.0005 0.0 0.0004 0.0005 0.0
-0.0004 0.0008 -0.0003 0.0013 0.0004 0.0014 0.0004 0.0002 -0.0001 0.0002 0.0001 -0.0002 -0.0002
-0.0008 -0.0015 0.0021 -0.0004 -0.0004 0.0012 0.0006 … -0.0011 0.0002 -0.0007 0.0001 0.0005 -0.0001
0.0008 -0.0004 0.0009 0.0019 -0.0022 -0.0014 0.0013 0.0003 -0.0001 -0.0003 -0.0 0.0003 0.0002
0.0005 0.0006 -0.0022 -0.0003 0.0 -0.0022 0.0016 -0.001 0.0002 -0.0003 0.0003 -0.0005 0.0005
0.001 -0.0005 -0.0007 0.0025 -0.0019 0.0005 0.0015 0.0002 -0.0005 0.0 -0.0005 0.0002 -0.0002
-0.0003 0.0002 -0.0021 -0.0008 -0.0012 0.0003 0.0003 -0.0001 -0.0 0.0001 0.0 0.0 0.0
Component Statistics:
Eigenvalues: 7975.439 4040.6348 2092.866 717.8178 468.5496 247.5029 202.2003 176.8452 87.3302 ... 0.3602
Cumulative Eigenvalues: 0.4779 0.72 0.8454 0.8884 0.9165 0.9313 0.9434 0.954 0.9593 ... 1.0
Explained SNR: 7974.4385 4039.6353 2091.8635 716.8176 467.55 246.5022 201.1985 175.845 86.3301 ... -0.6399
Cumulative SNR: 0.4839 0.729 0.856 0.8995 0.9278 0.9428 0.955 0.9657 0.9709 ... 0.9985
julia> transformed = forward_mnf(mnf, desis, 12);
julia> size(transformed)
(1131, 1120, 12)
julia> recovered = inverse_mnf(mnf, transformed);
julia> size(recovered)
(1131, 1120, 235)
RemoteSensingToolbox.cumulative_eigenvalues
— Methodcumulative_eigenvalues(x::MNF)
Return the cumulative eigenvalues associated with each principal components of the fitted MNF transform.
RemoteSensingToolbox.cumulative_snr
— Methodcumulative_snr(x::MNF)
Return the cumulative SNR associated with each principal components of the fitted MNF transform.
RemoteSensingToolbox.data_cov
— Methoddata_cov(x::MNF)
Return the data covariance matrix for the fitted MNF transform.
RemoteSensingToolbox.eigenvalues
— Methodeigenvalues(x::MNF)
Return the eigenvalues associated with each principal components of the fitted MNF transform.
RemoteSensingToolbox.estimate_noise
— Methodestimate_noise(raster::Union{<:AbstractRaster, <:AbstractRasterStack}; smooth=false)
Estimate the noise covariance matrix for the given raster.
Uses the Minimum/Maximum Autocorrelation Factor proposed by Switzer and Green.
For best results, the provided raster should be spectrally homoegenous (e.g., an open field or body of water).
Parameters
raster
: AnAbstractRaster
orAbstractRasterStack
.smooth
: Numerical stability requires that no bands have a variance of zero. A smoothing term can be applied to ensure that this is the case.
Example
using RemoteSensingToolbox, Rasters
# Load Data
src = DESIS("DESIS-HSI-L2A-DT0485529167_001-20220712T223540-V0220")
desis = decode(DESIS, Raster(src, :Bands))
# Extract Homogenous Region of Interest
roi = desis[X(1019:1040), Y(550:590)]
# Estimate Noise
ncm = estimate_noise(roi, smooth=true)
RemoteSensingToolbox.fit_mnf
— Methodfit_mnf(raster::Union{<:AbstractRasterStack, <:AbstractRaster}, noise_covariance::Matrix; stats_fraction=1.0)
fit_mnf(signatures::Matrix, noise_covariance::Matrix; stats_fraction=1.0)
Fit a Minimum Noise Fraction (MNF) transformation to the given AbstractRasterStack
or AbstractRaster
.
Parameters
raster
: TheAbstractRaster
orAbstractRasterStack
on which to fit the MNF transformation.signatures
: An n x b matrix of spectral signatures where n is the number of signatures and b is the number of bands.noise_sample
: A homogenous (spectrally uniform) region extracted fromraster
for calculating the noise covariance matrix.smooth
: The MNF transform cannot be computed if any band innoise_sample
has zero variance. To correct this, you may wish to introduce a small smoothing term (true by default).
RemoteSensingToolbox.forward_mnf
— Methodforward_mnf(transformation::MNF, raster, components::Int)
forward_mnf(transformation::MNF, signatures::Matrix, components::Int)
Perform a forward Minimum Noise Fraction (MNF) rotation on the given raster or signatures, retaining only the specified number of components.
Parameters
transformation
: A previously fitted MNF transformation.raster
: TheAbstractRaster
orAbstractRasterStack
on which to perform the MNF transformation.signatures
: An n x b matrix of spectral signatures where n is the number of signatures and b is the number of bands.components
: The number of bands to retain in the transformed image. All band numbers exceeding this value will be discarded.
RemoteSensingToolbox.inverse_mnf
— Methodinverse_mnf(transformation::MNF, raster::AbstractRaster)
inverse_mnf(transformation::MNF, signatures::Matrix)
Perform an inverse Minimum Noise Fraction (MNF) transformation to recover the original image or signatures.
Parameters
transformation
: A previously fitted MNF transformation.raster
: AnAbstractRaster
representing a previously transformed image. The number of bands should be less than or equal to that of the original image.signatures
: An n x p matrix of transformed signatures where n is the number of signatures and p is the number of retained components.
RemoteSensingToolbox.noise_cov
— Methodnoise_cov(x::MNF)
Return the noise covariance matrix for the fitted MNF transform.
RemoteSensingToolbox.projection
— Methodprojection(x::MNF)
Return the projection matrix for the fitted MNF transform.
RemoteSensingToolbox.snr
— Methodsnr(x::MNF)
Return the estimated SNR associated with each principal components of the fitted MNF transform.
Utilities
RemoteSensingToolbox
provides several utility functions for working with remotely sensed data.
RemoteSensingToolbox.apply_masks!
— Methodapply_masks!(raster, masks...)
Similar to Rasters.mask!
, but with the following differences:
- Removes non-missing mask values instead of missing values. This is useful when working with cloud or shadow masks.
- Accepts multiple masks, which are applied in sequence.
Parameters
raster
: TheAbstractRaster
orAbstractRasterStack
to be masked.masks
: One or more masks to apply to the given raster.
RemoteSensingToolbox.apply_masks
— Methodapply_masks(raster, masks...)
Similar to Rasters.mask
, but with the following differences:
- Removes non-missing mask values instead of missing values. This is useful when working with cloud or shadow masks.
- Accepts multiple masks, which are applied in sequence.
Parameters
raster
: TheAbstractRaster
orAbstractRasterStack
to be masked.masks
: One or more masks to apply to the given raster.
RemoteSensingToolbox.from_table
— Methodfrom_table(table, val_col, dim_col, dims; missingval=0)
Read a table into a Raster
or RasterStack
.
Parameters
table
: Any type that implements theTables.jl
interface.val_col
: The column(s) containing the raster values. Must be either a Symbol or Tuple of Symbols.dim_col
: The column(s) containing the coordinates of each value. Should be a Symbol if the coordinates
are stored as Tuples under a single column. It can also be a Tuple of Symbols if coordinates are stored as scalar values across multuple columns.
dims
: The dimensions of the output raster.missingval
: The value used to denote missing data.
Returns
Returns a single Raster
with the same dimensions as dims
when val_col
is a Symbol
. Otherwise, returns a RasterStack
in which each layer corresponds to a single value column.
RemoteSensingToolbox.has_bands
— Methodhas_bands(raster)
Returns true
if the provided raster or stack has a band dimension.
RemoteSensingToolbox.mask_nan!
— Methodmask_nan!(raster)
Replace all NaN
values with missingval(raster)
in-place.
Parameters
raster
: TheAbstractRaster
orAbstractRasterStack
from which we want to dropNaN
values.
RemoteSensingToolbox.nbands
— Methodnbands(raster)
Returns the number of spectral bands in the given AbstractRaster
or AbstractRasterStack
.
RemoteSensingToolbox.sample
— Methodsample(raster, [sink]; fraction=0.1)
Randomly sample a percentage of non-missing values from the provided raster.
Parameters
raster
: TheAbstractRaster
orAbstractRasterStack
from which to sample.sink
: ATables.jl
materializer (default =Tables.columntable
).fraction
: The fraction of values to sample (default = 10%).
Example
julia> src = Landsat8("LC08_L2SP_043024_20200802_20200914_02_T1/");
julia> rs = RasterStack(src, [:blue, :green, :red, :nir]);
julia> sample(rs, DataFrame)
4054017×4 DataFrame
Row │ B2 B3 B4 B5
│ UInt16 UInt16 UInt16 UInt16
─────────┼────────────────────────────────
1 │ 7841 9082 8174 24372
2 │ 8117 8977 8372 15577
3 │ 26460 26601 26579 27023
⋮ │ ⋮ ⋮ ⋮ ⋮
4054015 │ 7735 8403 7929 18386
4054016 │ 6984 9559 10026 11986
4054017 │ 8036 8658 8328 13896
4054011 rows omitted
RemoteSensingToolbox.statistics
— Methodstatistics(raster; stats=:all, fraction=1.0)
Calculate summary statistics for the provided Raster
or RasterStack
.
Parameters
raster
: AnAbstractRaster
orAbstractRasterStack
.stats
: AVector
orTuple
containing any combination of the symbols :mean, :std, and :cov.fraction
: The fraction of pixels to sample when calculating statistics.
Returns
A named tuple containing each requested statistic.
RemoteSensingToolbox.table
— Functiontable(raster, [sink])
Convert the raster into a table compatible with Tables.jl
.
Replaces all missingvals
with missing
.
Parameters
raster
: TheAbstractRaster
orAbstractRasterStack
to read into a table.sink
: ATables.jl
materializer (default =Tables.columntable
).
Example
julia> src = Landsat8("LC08_L2SP_043024_20200802_20200914_02_T1/");
julia> rs = RasterStack(src, [:blue, :green, :red, :nir]);
julia> table(rs, DataFrame) |> dropmissing!
40540174×5 DataFrame
Row │ geometry B2 B3 B4 B5
│ Tuple… UInt16 UInt16 UInt16 UInt16
──────────┼───────────────────────────────────────────────────────
1 │ (544335.0, 5.84578e6) 8345 8798 8216 14454
2 │ (544335.0, 5.84576e6) 8064 8707 8106 15583
3 │ (544365.0, 5.84576e6) 8247 8858 8135 17552
⋮ │ ⋮ ⋮ ⋮ ⋮ ⋮
40540172 │ (676365.0, 5.60968e6) 8863 9867 10164 16688
40540173 │ (676395.0, 5.60968e6) 8823 9684 10050 16210
40540174 │ (676425.0, 5.60968e6) 8934 9898 10324 16947
40540168 rows omitted