diff --git a/.gitignore b/.gitignore index 8077cf54..e0617f1e 100644 --- a/.gitignore +++ b/.gitignore @@ -7,3 +7,5 @@ Manifest.toml run-dev.sh .vscode/ benchmark/ +docs/build +docs/make-local.jl diff --git a/Project.toml b/Project.toml index 7a377aaf..5c78d416 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Omniscape" uuid = "a38d70fc-99f5-11e9-1e3d-cbca093024c3" authors = ["Vincent A. Landau "] -version = "0.5.4" +version = "0.5.5" [deps] ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3" @@ -16,7 +16,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] ArchGDAL = "0.4, 0.5" -Circuitscape = "5.9" +Circuitscape = "5.9.1" ProgressMeter = "1.3" StatsBase = "0.33" julia = "^1.5.4" diff --git a/docker/Dockerfile b/docker/Dockerfile index af547ae9..5e62778e 100644 --- a/docker/Dockerfile +++ b/docker/Dockerfile @@ -6,9 +6,9 @@ RUN apt-get update && \ curl # Install julia -RUN curl -L https://julialang-s3.julialang.org/bin/linux/x64/1.6/julia-1.6.2-linux-x86_64.tar.gz > /tmp/julia.tar.gz +RUN curl -L https://julialang-s3.julialang.org/bin/linux/x64/1.6/julia-1.6.3-linux-x86_64.tar.gz > /tmp/julia.tar.gz RUN tar -zxvf /tmp/julia.tar.gz -RUN ln -s /julia-1.6.2/bin/julia /usr/bin/julia +RUN ln -s /julia-1.6.3/bin/julia /usr/bin/julia # Install Omniscape RUN julia -e 'using Pkg; Pkg.add(["Omniscape", "Test", "BenchmarkTools"])' diff --git a/docs/src/apidocs.md b/docs/src/apidocs.md index 0fbc6bc1..49c27958 100644 --- a/docs/src/apidocs.md +++ b/docs/src/apidocs.md @@ -1,4 +1,8 @@ # API Documentation ```@docs run_omniscape + +missingarray + +missingarray_to_array ``` diff --git a/docs/src/examples.md b/docs/src/examples.md index 3691f964..36c752bc 100644 --- a/docs/src/examples.md +++ b/docs/src/examples.md @@ -71,7 +71,7 @@ reclass_table = [ ] ``` -Next, we define the configuration options for this model run. See the [Arguments](@ref) section in the [User Guide](@ref) for more information about each of the configuration options. +Next, we define the configuration options for this model run. See the [Settings and Options](@ref) section in the [User Guide](@ref) for more information about each of the configuration options. ```@example mdforest # Specify the configuration options diff --git a/src/Omniscape.jl b/src/Omniscape.jl index 7fba69fe..d7d3ba26 100644 --- a/src/Omniscape.jl +++ b/src/Omniscape.jl @@ -16,6 +16,6 @@ include("consts.jl") include("utils.jl") include("main.jl") include("errors_warnings.jl") -export run_omniscape +export run_omniscape, MissingArray, missingarray_to_array, missingarray end \ No newline at end of file diff --git a/src/config.jl b/src/config.jl index 1e8fb679..4e77b436 100644 --- a/src/config.jl +++ b/src/config.jl @@ -81,13 +81,15 @@ function init_csdict(cfg) a["ground_file_is_resistances"] = "True" a["use_direct_grounds"] = "False" - a["output_file"] = "temp" - a["write_cum_cur_map_only"] = "False" - + a["write_cum_cur_map_only"] = "False" a["scenario"] = "Advanced" - a["suppress_messages"] = cfg["suppress_cs_messages"] + a["connect_four_neighbors_only"] = cfg["connect_four_neighbors_only"] + a["cholmod_batch_size"] = "1000" + + check_solver!(cfg) + a["solver"] = cfg["solver"] a end diff --git a/src/main.jl b/src/main.jl index 65401304..87dbeb9e 100644 --- a/src/main.jl +++ b/src/main.jl @@ -7,21 +7,23 @@ In-memory method: run_omniscape( cfg::Dict{String, String} - resistance::Array{Union{Float64, Missing}, 2}; - reclass_table = Array{Union{Float64, Missing}, 2}(undef, 1, 2), + resistance::MissingArray{Number, 2}; + reclass_table = MissingArray{Float64, 2}(undef, 1, 2), source_strength = source_from_resistance(resistance, cfg, reclass_table), - condition1 = Array{Union{Float64, Missing}, 2}(undef, 1, 1), - condition2 = Array{Union{Float64, Missing}, 2}(undef, 1, 1), - condition1_future = Array{Union{Float64, Missing}, 2}(undef, 1, 1), - condition2_future = Array{Union{Float64, Missing}, 2}(undef, 1, 1), + condition1 = MissingArray{Float64, 2}(undef, 1, 1), + condition2 = MissingArray{Float64, 2}(undef, 1, 1), + condition1_future = MissingArray{Float64, 2}(undef, 1, 1), + condition2_future = MissingArray{Float64, 2}(undef, 1, 1), wkt = "", geotransform = [0.0, 1.0, 0.0, 0.0, 0.0, -1.0], write_outputs = false ) Compute omnidirectional current flow. All array inputs for the in-memory method -should be of type `Array{Union{T, Missing}, 2} where T <: Number`, with -`missing` used for NoData pixels. +should be of type `MissingArray{T, 2}`, which is just an alias for +`Array{Union{Missing, T}, 2}`. `missing` entries correspond to NoData pixels. +If you are starting with an array with a numeric NoData value, use +[`missingarray`](@ref) to convert it to a `MissingArray`. # Parameters **`path`**: The path to an INI file containing run parameters. See the @@ -92,13 +94,13 @@ key. """ function run_omniscape( cfg::Dict{String, String}, - resistance::Array{Union{T, Missing}, 2} where T <: Number; - reclass_table::Array{Union{T, Missing}, 2} where T <: Number = Array{Union{Float64, Missing}, 2}(undef, 1, 2), - source_strength::Array{Union{T, Missing}, 2} where T <: Number = source_from_resistance(resistance, cfg, reclass_table), - condition1::Array{Union{T, Missing}, 2} where T <: Number = Array{Union{Float64, Missing}, 2}(undef, 1, 1), - condition2::Array{Union{T, Missing}, 2} where T <: Number = Array{Union{Float64, Missing}, 2}(undef, 1, 1), - condition1_future::Array{Union{T, Missing}, 2} where T <: Number = condition1, - condition2_future::Array{Union{T, Missing}, 2} where T <: Number = condition2, + resistance::MissingArray{T, 2} where T <: Number; + reclass_table::MissingArray{T, 2} where T <: Number = MissingArray{Float64, 2}(undef, 1, 2), + source_strength::MissingArray{T, 2} where T <: Number = source_from_resistance(resistance, cfg, reclass_table), + condition1::MissingArray{T, 2} where T <: Number = MissingArray{Float64, 2}(undef, 1, 1), + condition2::MissingArray{T, 2} where T <: Number = MissingArray{Float64, 2}(undef, 1, 1), + condition1_future::MissingArray{T, 2} where T <: Number = condition1, + condition2_future::MissingArray{T, 2} where T <: Number = condition2, wkt::String = "", geotransform::Array{Float64, 1} = [0., 1., 0., 0., 0., -1.0], write_outputs::Bool = false) @@ -135,13 +137,9 @@ function run_omniscape( os_flags = get_omniscape_flags(cfg) # other - source_threshold = parse(Float64, cfg["source_threshold"]) project_name = cfg["project_name"] file_format = os_flags.write_as_tif ? "tif" : "asc" - check_solver!(cfg) - solver = cfg["solver"] - ## Set number of BLAS threads to 1 when parallel processing if os_flags.parallelize && nthreads() != 1 BLAS.set_num_threads(1) @@ -165,11 +163,22 @@ function run_omniscape( parse(Float64, cfg["condition2_lower"]), parse(Float64, cfg["condition2_upper"])) - condition_layers = ConditionLayers(condition1, condition1_future, condition2, condition2_future) + # a bit of a hack, but ensures compare_to_future can work in both methods + # for run_omniscape + if compare_to_future == "1" + condition2_future = condition2 + elseif compare_to_future == "2" + condition1_future = condition1 + elseif compare_to_future == "none" + condition1_future = condition1 + condition2_future = condition2 + end + + condition_layers = ConditionLayers(condition1, condition1_future, + condition2, condition2_future) ## Setup Circuitscape configuration cs_cfg_dict = init_csdict(cfg) - cs_cfg_dict["solver"] = solver cs_cfg = Circuitscape.init_config() Circuitscape.update!(cs_cfg, cs_cfg_dict) @@ -184,18 +193,11 @@ function run_omniscape( # Permute targets randomly (to get better estimate of ProgressMeter ETA) targets = targets[randperm(n_targets), :] - ## Define parameters for cs - # Get flags - o = Circuitscape.OutputFlags(false, false, - false, false, - false, false, - false, false) - precision_name = precision == Float64 ? "double" : "single" ## Add parallel workers if os_flags.parallelize @info("Starting up Omniscape with $(n_threads) workers and $(precision_name) precision") - @info("Using Circuitscape with the $(uppercase(solver)) solver...") + @info("Using Circuitscape with the $(cs_cfg["solver"]) solver...") cum_currmap = fill(convert(precision, 0.), int_arguments["nrows"], @@ -213,7 +215,7 @@ function run_omniscape( end else @info("Starting up Omniscape with 1 worker and $(precision_name) precision") - @info("Using Circuitscape with the $(uppercase(solver)) solver...") + @info("Using Circuitscape with the $(cs_cfg["solver"]) solver...") cum_currmap = fill(convert(precision, 0.), int_arguments["nrows"], int_arguments["ncols"], @@ -229,18 +231,11 @@ function run_omniscape( end end - cs_flags = Circuitscape.RasterFlags(true, false, true, - false, false, - false, Symbol("rmvsrc"), - cfg["connect_four_neighbors_only"] in TRUELIST, - false, solver, o) - if os_flags.correct_artifacts && !(int_arguments["block_size"] == 1) @info("Calculating block artifact correction array...") correction_array = calc_correction(int_arguments, os_flags, cs_cfg, - cs_flags, condition_layers, conditions, precision) @@ -272,7 +267,6 @@ function run_omniscape( resistance, os_flags, cs_cfg, - cs_flags, condition_layers, conditions, correction_array, @@ -291,18 +285,17 @@ function run_omniscape( target = Target(Int64(targets[i, 1]), Int64(targets[i, 2]), float(targets[i, 3])) try solve_target!(target, - int_arguments, - source_strength, - resistance, - os_flags, - cs_cfg, - cs_flags, - condition_layers, - conditions, - correction_array, - cum_currmap, - fp_cum_currmap, - precision) + int_arguments, + source_strength, + resistance, + os_flags, + cs_cfg, + condition_layers, + conditions, + correction_array, + cum_currmap, + fp_cum_currmap, + precision) catch error println("Omniscape failed on the moving window centered on row $(target.y_coord) column $(target.x_coord)") throw(error) @@ -361,9 +354,8 @@ function run_omniscape( write_cfg(cfg_user, project_name) if os_flags.reclassify && os_flags.write_reclassified_resistance - resistance[ismissing.(resistance)] .= -9999 write_raster("$(project_name)/classified_resistance", - convert(Array{precision, 2}, resistance), + missingarray_to_array(resistance, -9999), wkt, geotransform, file_format) @@ -403,17 +395,17 @@ function run_omniscape( end ## Return outputs, depending on user options # convert arrays, replace -9999's with missing - cum_currmap = convert_and_fill_missing(cum_currmap, precision) + cum_currmap = missingarray(cum_currmap, precision, -9999) if os_flags.calc_normalized_current && !os_flags.calc_flow_potential - normalized_cum_currmap = convert_and_fill_missing(normalized_cum_currmap, precision) + normalized_cum_currmap = missingarray(normalized_cum_currmap, precision, -9999) return cum_currmap, normalized_cum_currmap elseif !os_flags.calc_normalized_current && os_flags.calc_flow_potential - fp_cum_currmap = convert_and_fill_missing(fp_cum_currmap, precision) + fp_cum_currmap = missingarray(fp_cum_currmap, precision, -9999) return cum_currmap, fp_cum_currmap elseif os_flags.calc_normalized_current && os_flags.calc_flow_potential - fp_cum_currmap = convert_and_fill_missing(fp_cum_currmap, precision) - normalized_cum_currmap = convert_and_fill_missing(normalized_cum_currmap, precision) + fp_cum_currmap = missingarray(fp_cum_currmap, precision, -9999) + normalized_cum_currmap = missingarray(normalized_cum_currmap, precision, -9999) return cum_currmap, fp_cum_currmap, normalized_cum_currmap else return cum_currmap @@ -448,7 +440,7 @@ function run_omniscape(path::String) if os_flags.reclassify reclass_table = read_reclass_table("$(cfg["reclass_table"])", precision) else - reclass_table = Array{Union{precision, Missing}, 2}(undef, 1, 2) + reclass_table = MissingArray{precision, 2}(undef, 1, 2) end ## Load source strengths @@ -498,7 +490,7 @@ function run_omniscape(path::String) # get rid of unneeded raster to save memory condition1_future_raster = nothing else - condition1_future = condition1 + condition1_future = MissingArray{precision, 2}(undef, 1, 1) end if n_conditions == 2 @@ -525,16 +517,16 @@ function run_omniscape(path::String) # get rid of unneeded raster to save memory condition2_future_raster = nothing else - condition2_future = condition2 + condition2_future = MissingArray{precision, 2}(undef, 1, 1) end else - condition2 = Array{Union{Missing, precision}, 2}(undef, 1, 1) + condition2 = MissingArray{precision, 2}(undef, 1, 1) condition2_future = condition2 end else - condition1 = Array{Union{Missing, precision}, 2}(undef, 1, 1) - condition2 = Array{Union{Missing, precision}, 2}(undef, 1, 1) + condition1 = MissingArray{precision, 2}(undef, 1, 1) + condition2 = MissingArray{precision, 2}(undef, 1, 1) condition1_future = condition1 condition2_future = condition2 end diff --git a/src/structs.jl b/src/structs.jl index a49d2417..aeb22270 100644 --- a/src/structs.jl +++ b/src/structs.jl @@ -1,3 +1,8 @@ +import Base: size, getindex, setindex!, length, sum, eltype, view +import Base.Broadcast: broadcasted, broadcast +import StatsBase: mode +import Statistics: median + struct OmniscapeFlags calc_flow_potential::Bool calc_normalized_current::Bool @@ -30,10 +35,12 @@ struct Conditions condition2_upper::Number end -struct ConditionLayers - condition1_present::Array{Union{Missing, Number}, 2} # where T <: Number - condition1_future::Array{Union{Missing, Number}, 2} # where U <: Number - condition2_present::Array{Union{Missing, Number}, 2} # where V <: Number - condition2_future::Array{Union{Missing, Number}, 2} # where W <: Number -end +const MissingArray{T, N} = Array{Union{Missing, T}, N} + +struct ConditionLayers{T, N} + condition1_present::MissingArray{T, N} + condition1_future::MissingArray{T, N} + condition2_present::MissingArray{T, N} + condition2_future::MissingArray{T, N} +end diff --git a/src/utils.jl b/src/utils.jl index a52509b4..8c7d495c 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1,3 +1,5 @@ +import Circuitscape: compute_omniscape_current + function clip( A::Array{Union{Missing, T}, 2} where T <: Number; x::Int64, @@ -28,7 +30,6 @@ function clip( A_sub end - function get_targets( source_array::Array{Union{T, Missing}, 2} where T <: Number, arguments::Dict{String, Int64}, @@ -85,7 +86,7 @@ end # x and y defined by targets object. Ultimately the for loop will be done by # iterating through rows of targets object function get_source( - source_array::Array{Union{Missing, T}, 2} where T <: Number, + source_array::MissingArray{T, 2} where T <: Number, arguments::Dict{String, Int64}, os_flags::OmniscapeFlags, condition_layers::ConditionLayers, @@ -186,7 +187,7 @@ function get_source( source_subset end -function source_target_match!(source_subset::Array{Union{T, Missing}, 2} where T <: Number, +function source_target_match!(source_subset::MissingArray{T, 2} where T <: Number, n_conditions::Int64, condition_layers::ConditionLayers, conditions::Conditions, @@ -268,7 +269,7 @@ function get_ground(arguments::Dict{String, Int64}, end function get_conductance( - resistance::Array{Union{T, Missing}, 2} where T <: Number, + resistance::MissingArray{T, 2} where T <: Number, arguments::Dict{String, Int64}, target::Target, os_flags::OmniscapeFlags @@ -288,126 +289,13 @@ function get_conductance( convert(typeof(resistance_clipped), conductance) end - -function calculate_current( - conductance::Array{Union{T, Missing}, 2} where T <: Number, - source::Array{Union{T, Missing}, 2} where T <: Number, - ground::Array{T, 2} where T <: Number, - cs_flags::Circuitscape.RasterFlags, - cs_cfg::Dict{String, String}, - T::DataType - ) - V = Int64 - - # Replace missings with -9999, then convert to Array{T, 2} - # prior to circuitscape - conductance[ismissing.(conductance)] .= -9999 - source[ismissing.(source)] .= -9999 - conductance = convert(Array{T, 2}, conductance) - source = convert(Array{T, 2}, source) - - # get raster data - cellmap = conductance - polymap = Matrix{V}(undef, 0, 0) - source_map = source - ground_map = ground - points_rc = (V[], V[], V[]) - strengths = Matrix{T}(undef, 0, 0) - - included_pairs = Circuitscape.IncludeExcludePairs(:undef, - V[], - Matrix{V}(undef,0,0)) - - # This is just to satisfy type requirements, most of it not used - hbmeta = Circuitscape.RasterMeta(size(cellmap)[2], - size(cellmap)[1], - 0., - 0., - 1., - -9999., - Array{Float64, 1}(undef, 1), - "") - - rasterdata = Circuitscape.RasterData(cellmap, - polymap, - source_map, - ground_map, - points_rc, - strengths, - included_pairs, - hbmeta) - - # Generate advanced data - data = Circuitscape.compute_advanced_data(rasterdata, cs_flags, cs_cfg) - - G = data.G - nodemap = data.nodemap - polymap = data.polymap - hbmeta = data.hbmeta - sources = data.sources - grounds = data.grounds - finitegrounds = data.finitegrounds - cc = data.cc - check_node = data.check_node - source_map = data.source_map # Need it for one to all mode - cellmap = data.cellmap - - f_local = Vector{eltype(G)}() - voltages = Vector{eltype(G)}() - outcurr = Circuitscape.alloc_map(hbmeta) - - for c in cc - if check_node != -1 && !(check_node in c) - continue - end - - # a_local = laplacian(G[c, c]) - a_local = G[c,c] - s_local = sources[c] - g_local = grounds[c] - - if sum(s_local) == 0 || sum(g_local) == 0 - continue - end - - if finitegrounds != [-9999.] - f_local = finitegrounds[c] - else - f_local = finitegrounds - end - - voltages = Circuitscape.multiple_solver(cs_cfg, - data.solver, - a_local, - s_local, - g_local, - f_local) - - local_nodemap = Circuitscape.construct_local_node_map(nodemap, - c, - polymap) - - Circuitscape.accum_currents!(outcurr, - voltages, - cs_cfg, - a_local, - voltages, - f_local, - local_nodemap, - hbmeta) - end - - outcurr -end - function solve_target!( target::Target, int_arguments::Dict{String, Int64}, - source_strength::Array{Union{Missing, T}, 2} where T <: Number, - resistance::Array{Union{Missing, T}, 2} where T <: Number, + source_strength::MissingArray{T, 2} where T <: Number, + resistance::MissingArray{T, 2} where T <: Number, os_flags::OmniscapeFlags, cs_cfg::Dict{String, String}, - cs_flags::Circuitscape.RasterFlags, condition_layers::ConditionLayers, conditions::Conditions, correction_array::Array{T, 2} where T <: Number, @@ -438,23 +326,22 @@ function solve_target!( grid_size = size(source) ## Run circuitscape - curr = calculate_current(conductance, - source, - ground, - cs_flags, - cs_cfg, - precision) + conductance = missingarray_to_array(conductance, -9999) + source = missingarray_to_array(source, -9999) + + curr = compute_omniscape_current(conductance, + source, + ground, + cs_cfg) ## If normalize = True, calculate null map and normalize if os_flags.compute_flow_potential - null_conductance = convert(Array{Union{precision, Missing}, 2}, fill(1, grid_size)) - - flow_potential = calculate_current(null_conductance, - source, - ground, - cs_flags, - cs_cfg, - precision) + null_conductance = convert(Array{precision, 2}, fill(1, grid_size)) + + flow_potential = compute_omniscape_current(null_conductance, + source, + ground, + cs_cfg) end if os_flags.correct_artifacts && !(int_arguments["block_size"] == 1) @@ -512,7 +399,6 @@ function calc_correction( arguments::Dict{String, Int64}, os_flags::OmniscapeFlags, cs_cfg::Dict{String, String}, - cs_flags::Circuitscape.RasterFlags, condition_layers::ConditionLayers, conditions::Conditions, precision::DataType @@ -522,18 +408,23 @@ function calc_correction( # are not adjusted by target weight, but stay the same according to their # original values. Something to keep in mind... - temp_source = convert(Array{Union{precision, Missing}, 2}, - fill(1.0, - arguments["radius"] * 2 + buffer * 2 + 1, - arguments["radius"] * 2 + buffer * 2 + 1)) + temp_source = convert( + Array{precision, 2}, + fill( + 1.0, + arguments["radius"] * 2 + buffer * 2 + 1, + arguments["radius"] * 2 + buffer * 2 + 1 + ) + ) + temp_source = missingarray(temp_source, precision, -9999) source_null = clip(temp_source, x = arguments["radius"] + buffer + 1, y = arguments["radius"] + buffer + 1, distance = arguments["radius"]) - # Append NoData (-9999) if buffer > 0 - if buffer > 0 + # Append NoData (missing) if buffer > 0 + if buffer > 0 column_dims = (size(source_null)[1], buffer) # Add columns source_null = hcat(fill(missing, column_dims), @@ -546,7 +437,8 @@ function calc_correction( source_null, fill(missing, row_dims)) end - n_sources = sum(source_null[(!).(ismissing.(source_null))]) + + n_sources = sum(skipmissing(source_null)) source_null[ismissing.(source_null)] .= 0.0 source_null[arguments["radius"] + arguments["buffer"] + 1, @@ -576,19 +468,21 @@ function calc_correction( arguments["radius"] + arguments["buffer"] + 1] = Inf - block_null_current = calculate_current(conductance, - source_blocked, - ground, - cs_flags, - cs_cfg, - precision) + # Convert inputs for Circuitscape current solve + conductance = missingarray_to_array(conductance, -9999) + source_blocked = missingarray_to_array(source_blocked, -9999) + source_null = missingarray_to_array(source_null, -9999) + + block_null_current = compute_omniscape_current(conductance, + source_blocked, + ground, + cs_cfg) + + null_current = compute_omniscape_current(conductance, + source_null, + ground, + cs_cfg) - null_current = calculate_current(conductance, - source_null, - ground, - cs_flags, - cs_cfg, - precision) null_current_total = fill(convert(precision, 0.), arguments["radius"] * 2 + arguments["buffer"] * 2 + arguments["block_size"], arguments["radius"] * 2 + arguments["buffer"] * 2 + arguments["block_size"]) @@ -633,9 +527,9 @@ function get_omniscape_flags(cfg::Dict{String, String}) end # Calculate the source layer using resistance surface and arguments from cfg -function source_from_resistance(resistance::Array{Union{T, Missing}, 2} where T <: Number, +function source_from_resistance(resistance::MissingArray{T, 2} where T <: Number, cfg::Dict{String, String}, - reclass_table::Array{Union{T, Missing}, 2} where T <: Number) + reclass_table::MissingArray{T, 2} where T <: Number) full_cfg = init_cfg() update_cfg!(full_cfg, cfg) r_cutoff = parse(Float64, full_cfg["r_cutoff"]) @@ -659,9 +553,8 @@ function source_from_resistance(resistance::Array{Union{T, Missing}, 2} where T source_strength end - -function reclassify_resistance!(resistance::Array{Union{T, Missing}, 2} where T <: Number, - reclass_table::Array{Union{T, Missing}, 2} where T <: Number) +function reclassify_resistance!(resistance::MissingArray{T, 2} where T <: Number, + reclass_table::MissingArray{T, 2} where T <: Number) resistance_old = deepcopy(resistance) for i in 1:(size(reclass_table)[1]) resistance[coalesce.(resistance_old .== reclass_table[i, 1], false)] .= reclass_table[i, 2] @@ -669,19 +562,62 @@ function reclassify_resistance!(resistance::Array{Union{T, Missing}, 2} where T resistance_old = nothing # remove from memory end -function convert_and_fill_missing(A::Array{T, 2} where T <: Number, - precision::DataType) - A = convert(Array{Union{precision, Missing}, 2}, A) - A[A .== -9999] .= missing - - A -end - -function arrays_equal(A::Array{Union{T, Missing}, 2} where T <: Number, - B::Array{Union{T, Missing}, 2} where T <: Number) +function arrays_equal(A::MissingArray{T, 2} where T <: Number, + B::MissingArray{T, 2} where T <: Number) # Check that non-missing entries are equal A[ismissing.(A)] .= -9999 B[ismissing.(B)] .= -9999 isapprox(Array{Float64, 2}(A), Array{Float64, 2}(B); rtol = 1e-6) end + +""" + missingarray(A::Array{T, N}, T::DataType, nodata::Number) + +This function converts an array to a `MissingArray` and replaces `nodata` +values with `missing` in the output. `MissingArray{T, N}` is an alias +for `Array{Union{T, Missing}, N}``. This function can be used to prepare +inputs for [`run_omniscape`](@ref). + +# Parameters + +**`A`**: The array to convert. + +**`T`**: The data type for the output (e.g. Float64 or Float32). + +**`nodata`**: The numeric value to be replaced by `missing` in the result. +""" +function missingarray( + A::Array{T, N} where T <: Union{Missing, Number} where N, + T::DataType, + nodata::Number + ) + output = convert(Array{Union{T, Missing}, ndims(A)}, copy(A)) + output[output .== nodata] .= missing + + return output +end + +""" + missing_array_to_array(A::MissingArray{T, N}, nodata::Number) + +This function converts an array of type `MissingArray` to a numeric array +and replaces `missing` entries with `nodata`. `MissingArray{T, N}` is an alias +for `Array{Union{T, Missing}, N}`. + +# Parameters + +**`A`**: The array to convert. + +**`nodata`**: The numeric value with which `missing` values will be replaced in +the result. +""" +function missingarray_to_array( + A::MissingArray{T, N} where T <: Number where N, + nodata::Number + ) + output = copy(A) + output[ismissing.(output)] .= nodata + + return convert(Array{typeof(output[1]), ndims(output)}, output) +end