Skip to content

Commit

Permalink
Issue #25 - Add code documentation
Browse files Browse the repository at this point in the history 
- Add docstrings and apply formatting to:
    - src/EndmemberLibrary.jl
    - src/Plotting.jl
    - src/Solvers.jl
    - src/Datasets.jl
    - src/SpectralUnmixing.jl

- Minor README.md updates

- Update example notebook with endmember library citation
  • Loading branch information
Vatsal A Jhalani committed Nov 4, 2024
1 parent 40c9862 commit 9bbef4d
Show file tree
Hide file tree
Showing 7 changed files with 1,008 additions and 251 deletions.
9 changes: 5 additions & 4 deletions README.md
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
<h1 align="center">
SpectralUnmixing
</h1>
A general, fast, flexible, and including spectral unmixing package. Oriented towards VSWIR imaging spectroscopy data but applicable for different sensor types. Includes options for different treatments of endmember library assemblages, including MESMA and bootstrapping (aka monte carlo) strategies.
A general, fast, flexible, and including spectral unmixing package. Oriented towards VSWIR imaging spectroscopy data but applicable for different sensor types. Includes options for different treatments of endmember library assemblages, including MESMA and bootstrapping (aka Monte Carlo) strategies.

## Installation
This package has not been registered yet, but will be soon. In the interim, after cloning and navigating into the repository, it can be installed from the Julia REPL, by running
Expand All @@ -12,17 +12,19 @@ export JULIA_PROJECT=${PWD}
```

## Using the script
Currently the package supports reading and writing ENVI raster data.

Basic:

```
julia unmix.jl REFLECTANCE_IMAGE ENDMEMBER_LIBRARY ENDMEMBER_COLUMN OUTPUT_BASE --mode sma
julia unmix.jl REFLECTANCE_IMAGE ENDMEMBER_LIBRARY ENDMEMBER_COLUMN OUTPUT_BASE --mode sma
```


Parallel implementation (with 10 cores):

```
julia -p 10 unmix.jl REFLECTANCE_IMAGE ENDMEMBER_LIBRARY ENDMEMBER_COLUMN OUTPUT_BASE --mode sma
julia -p 10 unmix.jl REFLECTANCE_IMAGE ENDMEMBER_LIBRARY ENDMEMBER_COLUMN OUTPUT_BASE --mode sma
```

Bootstrapping uncertainty:
Expand All @@ -42,4 +44,3 @@ Preset maximum number of endmembers used for unmixing:
```
julia -p 10 unmix.jl REFLECTANCE_IMAGE ENDMEMBER_LIBRARY ENDMEMBER_COLUMN OUTPUT_BASE --mode sma --n_mc 50 --normalization brightness --num_endmembers 10
```

2 changes: 1 addition & 1 deletion examples/simulation_demo.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,7 @@
"source": [
"# Spectral Library Preparation\n",
"\n",
"Start by loading up a spectral library, filtering it down to the desired classes, interpolating the spectra to the relevant wavelengths, and then removing undesired wavelength regions. By default the ignored regions correspond to atmospheric absorption regions. Here we split the library into the highest level classes: soil, photosynthetic vegetation (pv) and non-photosynthetic vegetation (npv). As the plots below show, all the endmember spectra within each class share general attributes while having small variations between them."
"Start by loading up a spectral library (Ochoa et al. [10.21232/QijXnpFt](https://doi.org/10.21232/QijXnpFt)), filtering it down to the desired classes, interpolating the spectra to the relevant wavelengths, and then removing undesired wavelength regions. By default the ignored regions correspond to atmospheric absorption regions. Here we split the library into the highest level classes: soil, photosynthetic vegetation (pv) and non-photosynthetic vegetation (npv). As the plots below show, all the endmember spectra within each class share general attributes while having small variations between them."
]
},
{
Expand Down
245 changes: 191 additions & 54 deletions src/Datasets.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,17 @@
using ArchGDAL
using GDAL

"""
set_band_names(filename::String, band_names::Vector)
Set the names of the bands in a raster dataset specified by the given filename.
# Arguments
- `filename::String`: Path to the raster file. The file must be in a format supported by
GDAL. The data file is updated in place.
- `band_names::Vector`: A vector of strings containing the new names for each band in the
raster dataset. The length of this vector must match the number of bands in the raster.
"""
function set_band_names(filename::String, band_names::Vector)
GC.gc()
local_ods = GDAL.gdalopen(filename, GDAL.GA_Update)
Expand All @@ -27,45 +37,104 @@ function set_band_names(filename::String, band_names::Vector)
local_ods = nothing
end

function initiate_output_datasets(output_files::Vector{String}, x_len::Int64, y_len::Int64, output_bands::Vector, reference_dataset::ArchGDAL.IDataset)
"""
initiate_output_datasets(output_files::Vector{String}, x_len::Int64, y_len::Int64,
output_bands::Vector, reference_dataset::ArchGDAL.IDataset)
Initialize multiple output raster datasets in ENVI format with specified dimensions and
band information using a reference dataset for projection and geotransformation.
# Arguments
- `output_files::Vector{String}`: File paths of the output datasets. Each entry corresponds
to a separate output raster.
- `x_len::Int64`: Width of the output datasets.
- `y_len::Int64`: Height of the output datasets.
- `output_bands::Vector`: Number of bands for each output dataset. The length of this
vector must match the length of `output_files`.
- `reference_dataset::ArchGDAL.IDataset`: Reference dataset object from which projection
and geotransformation information will be copied.
# Notes
- All bands in the output datasets are initialized with the no-data value of `-9999`.
"""
function initiate_output_datasets(output_files::Vector{String}, x_len::Int64, y_len::Int64,
output_bands::Vector, reference_dataset::ArchGDAL.IDataset)
GC.gc()
for _o in 1:length(output_files)
outDataset = ArchGDAL.create(output_files[_o], driver=ArchGDAL.getdriver("ENVI"), width=x_len,
height=y_len, nbands=output_bands[_o], dtype=Float64)
outDataset = ArchGDAL.create(
output_files[_o], driver=ArchGDAL.getdriver("ENVI"), width=x_len,
height=y_len, nbands=output_bands[_o], dtype=Float64
)
ArchGDAL.setproj!(outDataset, ArchGDAL.getproj(reference_dataset))
@info "Output Image Size (x,y,b): $x_len, $y_len, $output_bands. Creating output fractional cover dataset."
@info "Output Image Size (x,y,b): $x_len, $y_len, $output_bands.
Creating output fractional cover dataset."
try
ArchGDAL.setgeotransform!(outDataset, ArchGDAL.getgeotransform(reference_dataset))
ArchGDAL.setgeotransform!(
outDataset, ArchGDAL.getgeotransform(reference_dataset)
)
catch e
println("No geotransorm avaialble, proceeding without")
println("No geotransorm available, proceeding without")
end

for _b in 1:length(output_bands)
ArchGDAL.setnodatavalue!(ArchGDAL.getband(outDataset,_b), -9999)
ArchGDAL.setnodatavalue!(ArchGDAL.getband(outDataset, _b), -9999)
end

output = zeros(x_len, y_len, output_bands[_o]) .- 9999
ArchGDAL.write!(outDataset, output, [1:size(output)[end];], 0, 0, size(output)[1], size(output)[2])
ArchGDAL.write!(
outDataset, output, [1:size(output)[end];], 0, 0, size(output)[1],
size(output)[2]
)

outDataset = nothing

end
end

"""
write_results(output_files::Vector{String}, output_bands::Vector, x_len::Int64,
y_len::Int64, results, args)
Write processed results to specified output ENVI raster files, including primary outputs,
uncertainty estimates, and complete fractions, based on the provided parameters.
# Arguments
- `output_files::Vector{String}`: File paths of the output datasets to be created or
updated. Each entry corresponds to a separate output raster.
- `output_bands::Vector`: Number of bands for each output dataset. The length of this
vector must match the number of output files.
- `x_len::Int64`: Width of the output datasets.
- `y_len::Int64`: Height of the output datasets.
- `results`: A collection of results, where each result is expected to be a tuple or
similar structure. Specific elements contain the data to be written to the output datasets.
- `args`: An object or structure containing additional parameters that dictate the
writing behavior. Available options are:
- `args.n_mc`: If greater than 1, an uncertainty output is also written to the second
output file.
- `args.write_complete_fractions`: If set to 1, write complete fractions to the
subsequent output file.
# Notes
- The primary output is initialized with a no-data value of `-9999` and written to the
first output file specified in `output_files`.
- Each output array is permuted to match the desired dimensions before writing.
"""
function write_results(output_files::Vector{String}, output_bands::Vector, x_len::Int64,
y_len::Int64, results, args)

function write_results(output_files::Vector{String}, output_bands::Vector, x_len::Int64, y_len::Int64, results, args)
GC.gc()

@info "Write primary output"
output = zeros(y_len, x_len, output_bands[1]) .- 9999
for res in results
if isnothing(res[2]) == false
output[res[1],res[3], :] = res[2]
output[res[1], res[3], :] = res[2]
end
end
output = permutedims( output, (2,1,3))
outDataset = ArchGDAL.read(output_files[1], flags=1, alloweddrivers =["ENVI"])
ArchGDAL.write!(outDataset, output, [1:size(output)[end];], 0, 0, size(output)[1], size(output)[2])
output = permutedims(output, (2, 1, 3))
outDataset = ArchGDAL.read(output_files[1], flags=1, alloweddrivers=["ENVI"])
ArchGDAL.write!(
outDataset, output, [1:size(output)[end];], 0, 0, size(output)[1], size(output)[2]
)
outDataset = nothing

ods_idx = 2
Expand All @@ -74,13 +143,16 @@ function write_results(output_files::Vector{String}, output_bands::Vector, x_len
output = zeros(y_len, x_len, output_bands[ods_idx]) .- 9999
for res in results
if isnothing(res[4]) == false
output[res[1],res[3], :] = res[4]
output[res[1], res[3], :] = res[4]
end
end

output = permutedims( output, (2,1,3))
outDataset = ArchGDAL.read(output_files[ods_idx], flags=1, alloweddrivers =["ENVI"])
ArchGDAL.write!(outDataset, output, [1:size(output)[end];], 0, 0, size(output)[1], size(output)[2])
output = permutedims(output, (2, 1, 3))
outDataset = ArchGDAL.read(output_files[ods_idx], flags=1, alloweddrivers=["ENVI"])
ArchGDAL.write!(
outDataset, output, [1:size(output)[end];], 0, 0, size(output)[1],
size(output)[2]
)
outDataset = nothing
ods_idx += 1
end
Expand All @@ -90,29 +162,53 @@ function write_results(output_files::Vector{String}, output_bands::Vector, x_len
output = zeros(y_len, x_len, output_bands[ods_idx]) .- 9999
for res in results
if isnothing(res[5]) == false
output[res[1],res[3], :] = res[5]
output[res[1], res[3], :] = res[5]
end
end

output = permutedims( output, (2,1,3))
outDataset = ArchGDAL.read(output_files[ods_idx], flags=1, alloweddrivers =["ENVI"])
ArchGDAL.write!(outDataset, output, [1:size(output)[end];], 0, 0, size(output)[1], size(output)[2])
output = permutedims(output, (2, 1, 3))
outDataset = ArchGDAL.read(output_files[ods_idx], flags=1, alloweddrivers=["ENVI"])
ArchGDAL.write!(
outDataset, output, [1:size(output)[end];], 0, 0, size(output)[1],
size(output)[2]
)
outDataset = nothing
ods_idx += 1
end

end

function write_line_results(output_files::Vector{String}, results, n_mc::Int64, write_complete_fractions::Bool)
"""
write_line_results(output_files::Vector{String}, results, n_mc::Int64,
write_complete_fractions::Bool)
Write line-based results to specified output ENVI raster files, including primary results,
uncertainty estimates, and complete fractions based on provided parameters.
# Arguments
- `output_files::Vector{String}`: File paths of the output datasets to be created or
updated. Each entry corresponds to a separate output raster.
- `results`: A collection of results, where each entry is expected to be a tuple or
similar structure. Specific elements contain the data to be written to the output datasets.
- `n_mc::Int64`: If greater than one, uncertainty outputs will be written to the
subsequent output file.
- `write_complete_fractions::Bool`: Indicates whether to write complete fractions
to the subsequent output file.
"""
function write_line_results(output_files::Vector{String}, results, n_mc::Int64,
write_complete_fractions::Bool)

GC.gc()

if isnothing(results[2]) == false

output = zeros(size(results[3])[1], 1, size(results[2])[2]) .- 9999
output[results[3], 1, :] = results[2]

outDataset = ArchGDAL.read(output_files[1], flags=1, alloweddrivers =["ENVI"])
ArchGDAL.write!(outDataset, output, [1:size(output)[end];], 0, results[1]-1, size(output)[1], 1)
outDataset = ArchGDAL.read(output_files[1], flags=1, alloweddrivers=["ENVI"])
ArchGDAL.write!(
outDataset, output, [1:size(output)[end];], 0, results[1] - 1,
size(output)[1], 1
)
outDataset = nothing
end
ods_idx = 2
Expand All @@ -121,9 +217,14 @@ function write_line_results(output_files::Vector{String}, results, n_mc::Int64,
if isnothing(results[4]) == false
output = zeros(size(results[3])[1], 1, size(results[4])[2]) .- 9999
output[results[3], 1, :] = results[4]

outDataset = ArchGDAL.read(output_files[ods_idx], flags=1, alloweddrivers =["ENVI"])
ArchGDAL.write!(outDataset, output, [1:size(output)[end];], 0, results[1]-1, size(output)[1], 1)

outDataset = ArchGDAL.read(
output_files[ods_idx], flags=1, alloweddrivers=["ENVI"]
)
ArchGDAL.write!(
outDataset, output, [1:size(output)[end];], 0, results[1] - 1,
size(output)[1], 1
)
outDataset = nothing
end
ods_idx += 1
Expand All @@ -133,36 +234,72 @@ function write_line_results(output_files::Vector{String}, results, n_mc::Int64,
if isnothing(results[5]) == false
output = zeros(size(results[3])[1], 1, size(results[5])[2]) .- 9999
output[results[3], 1, :] = results[5]

outDataset = ArchGDAL.read(output_files[ods_idx], flags=1, alloweddrivers =["ENVI"])
ArchGDAL.write!(outDataset, output, [1:size(output)[end];], 0, results[1]-1, size(output)[1], 1)

outDataset = ArchGDAL.read(
output_files[ods_idx], flags=1, alloweddrivers=["ENVI"]
)
ArchGDAL.write!(
outDataset, output, [1:size(output)[end];], 0, results[1] - 1,
size(output)[1], 1
)
outDataset = nothing
end
ods_idx += 1
end

end

function load_line(reflectance_file::String, reflectance_uncertainty_file::String, line::Int64,
good_bands::Array{Bool}, refl_nodata::Float64)

reflectance_dataset = ArchGDAL.read(reflectance_file, alloweddrivers =["ENVI"])
img_dat = convert(Array{Float64},ArchGDAL.readraster(reflectance_file, alloweddrivers =["ENVI"])[:,line,:])
img_dat = img_dat[:, good_bands]
good_data = .!all(img_dat .== refl_nodata, dims=2)[:,1]
img_dat = img_dat[good_data,:]

if sum(good_data) >= 1
if reflectance_uncertainty_file != ""
unc_dat = convert(Array{Float64},ArchGDAL.readraster(reflectance_uncertainty_file, alloweddrivers =["ENVI"])[:,line,:])
unc_dat = unc_dat[:, good_bands]
unc_dat = unc_dat[good_data,:]
else
unc_dat = nothing
end
"""
load_line(reflectance_file::String, reflectance_uncertainty_file::String,
line::Int64, good_bands::Array{Bool}, refl_nodata::Float64)
Load a specific line (row) of reflectance data and its associated uncertainty from raster
files, filtering based on specified good bands and no-data values.
# Arguments
- `reflectance_file::String`: Path to the ENVI reflectance raster file.
- `reflectance_uncertainty_file::String`: Path to the uncertainty raster file corresponding
to the reflectance data. An empty string indicates no uncertainty data.
- `line::Int64`: Line (row) of the raster data to be loaded.
- `good_bands::Array{Bool}`: Boolean array indicating which bands of the data to keep and
process.
- `refl_nodata::Float64`: The no-data value used in the reflectance data.
# Returns
- A tuple containing:
- `img_dat::Array{Float64}`: The filtered reflectance data for the specified line and
good bands. If no valid data exists, this will be `nothing`.
- `unc_dat::Array{Float64}`: The uncertainty data, filtered similarly. If no uncertainty
data is loaded or valid, this will be `nothing`.
- `good_data`: A boolean array indicating which pixels contain valid reflectance data.
"""
function load_line(reflectance_file::String, reflectance_uncertainty_file::String,
line::Int64, good_bands::Array{Bool}, refl_nodata::Float64)

reflectance_dataset = ArchGDAL.read(reflectance_file, alloweddrivers=["ENVI"])
img_dat = convert(
Array{Float64},
ArchGDAL.readraster(reflectance_file, alloweddrivers=["ENVI"])[:, line, :]
)
img_dat = img_dat[:, good_bands]
good_data = .!all(img_dat .== refl_nodata, dims=2)[:, 1]
img_dat = img_dat[good_data, :]

if sum(good_data) >= 1
if reflectance_uncertainty_file != ""
unc_dat = convert(
Array{Float64},
ArchGDAL.readraster(
reflectance_uncertainty_file, alloweddrivers=["ENVI"]
)[:, line, :]
)
unc_dat = unc_dat[:, good_bands]
unc_dat = unc_dat[good_data, :]
else
return nothing, nothing, good_data
unc_dat = nothing
end
else
return nothing, nothing, good_data
end

return img_dat, unc_dat, good_data
return img_dat, unc_dat, good_data
end
Loading

0 comments on commit 9bbef4d

Please sign in to comment.