Use @lock for adding moving window results to cumulative current array

[vlandau]

This will result in additional overhead, but significant memory reductions because the cumulative current storage array won't need to have an additional dimension of length n = <number of threads>, it can just be stored in a matrix. I imagine this should be done in addition to #79, where we use a hierarchical parallel processing framework with chunks of the landscape parallelized across processes, and individual moving window solves parallelized via multi-threading.

To start, I'll do some basic benchmarking at different numbers of threads with the @lock method to see what level of overhead we're talking here. Might be minimal, but I imagine it will increase as the number of threads increases because we'll have a "longer line" of processes waiting to write to the array.

[vlandau]

Okay -- this was easy to implement. I've got two methods side y side to allow benchmarking -- still need to do a formal benchmark to see how this scales with increased numbers of threads. Here's the relevant commit.

Here's a script to try out both methods. With 6 threads, they're super similar in compute time, but the lock method has the benefit of much lower memory usage.

# using Pkg; Pkg.add(["GeoData", "Plots"])
using Omniscape

url_base = "https://raw.githubusercontent.com/Circuitscape/datasets/main/"
# Download the NLCD tile used to create the resistance surface and load it
download(string(url_base, "data/nlcd_2016_frederick_md.tif"),
         "nlcd_2016_frederick_md.tif")

# # Plot the landcover data
# values = [11, 21, 22, 23, 24, 31, 41, 42, 43, 52, 71, 81, 82, 90, 95]
# palette = ["#476BA0", "#DDC9C9", "#D89382", "#ED0000", "#AA0000",
#            "#b2b2b2", "#68AA63", "#1C6330", "#B5C98E", "#CCBA7C",
#            "#E2E2C1", "#DBD83D", "#AA7028", "#BAD8EA", "#70A3BA"]

# plot(GeoArray(GDALarray("nlcd_2016_frederick_md.tif")),
#      title = "Land Cover Type", xlabel = "Easting", ylabel = "Northing",
#      seriescolor = cgrad(palette, (values .- 12) ./ 84, categorical = true),
#      size = (700, 640))

land_cover, wkt, transform = Omniscape.read_raster("nlcd_2016_frederick_md.tif", Float64)


# Create the reclassification table used to translate land cover into resistance
reclass_table = [
    11.	100; # Water
    21	500; # Developed, open space
    22	1000; # Developed, low intensity
    23	missing; # Developed, medium intensity
    24	missing; # Developed, high intensity
    31	100; # Barren land
    41	1; # Deciduous forest
    42	1; # Evergreen forest
    43	1; # Mixed forest
    52	20; # Shrub/scrub
    71	30; # Grassland/herbaceous
    81	200; # Pasture/hay
    82	300; # Cultivated crops
    90	20; # Woody wetlands
    95	30; # Emergent herbaceous wetlands
]


config = Dict{String, String}(
    "radius" => "100",
    "block_size" => "21",
    "project_name" => "md_nlcd_omniscape_output",
    "source_from_resistance" => "true",
    "r_cutoff" => "1", # Only forest pixels should be sources
    "reclassify_resistance" => "true",
    "calc_normalized_current" => "true",
    "calc_flow_potential" => "true"
)

currmap, flow_pot, norm_current = run_omniscape_lock(config,
                                                land_cover,
                                                reclass_table = reclass_table,
                                                wkt = wkt,
                                                geotransform = transform,
                                                write_outputs = true)

currmap, flow_pot, norm_current = run_omniscape(config,
                                                land_cover,
                                                reclass_table = reclass_table,
                                                wkt = wkt,
                                                geotransform = transform,
                                                write_outputs = true)