Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Issue #25 - Add code documentation #26

Merged
merged 1 commit into from
Dec 19, 2024
Merged

Issue #25 - Add code documentation #26

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 5 additions & 4 deletions README.md
View file
Open in desktop
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
```

245 changes: 191 additions & 54 deletions src/Datasets.jl
View file
Open in desktop
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.
"""
pgbrodrick marked this conversation as resolved.
Show resolved Hide resolved
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
Loading