feat: data initialisation, criteria and assessment parameters refactor

src/ReefGuideAPI.jl

@@ -1,93 +1,50 @@
 module ReefGuideAPI
 
+# System imports 
 using Base.Threads
-using
-    Glob,
-    TOML
 
-using Serialization
+# File IO
+using Glob, TOML, Serialization
 
-using DataFrames
-using OrderedCollections
-using Memoization
-using SparseArrays
+# Geospatial
+using ArchGDAL, GeoParquet, Rasters
 
-using FLoops, ThreadsX
-
-import GeoDataFrames as GDF
-using
-    ArchGDAL,
-    GeoParquet,
-    Rasters
-
-using
-    HTTP,
-    Oxygen
-
-include("job_worker/Worker.jl")
+# Server
+using HTTP, Oxygen
 
-include("Middleware.jl")
-include("admin.jl")
-include("file_io.jl")
-include("server_cache.jl")
+# Collections
+using DataFrames, OrderedCollections, SparseArrays
 
-# TODO Remove these due to deprecation
-include("job_management/JobInterface.jl")
-include("job_management/DiskService.jl")
-
-include("criteria_assessment/criteria.jl")
-include("criteria_assessment/query_thresholds.jl")
-include("criteria_assessment/regional_assessment.jl")
-
-include("site_assessment/common_functions.jl")
-include("site_assessment/best_fit_polygons.jl")
+# Multithreading
+using FLoops, ThreadsX
 
-function get_regions()
-    # TODO: Comes from config?
-    regions = String[
-        "Townsville-Whitsunday",
-        "Cairns-Cooktown",
-        "Mackay-Capricorn",
-        "FarNorthern"
-    ]
+# Utilities and helpers for assessments
+include("utility/utility.jl")
 
-    return regions
-end
+# Assessment logic
+include("assessment_methods/assessment_methods.jl")
 
-function get_auth_router(config::Dict)
-    # Setup auth middleware - depends on config.toml - can return identity func
-    auth = setup_jwt_middleware(config)
-    return router(""; middleware=[auth])
-end
+# Worker system
+include("job_worker/job_worker.jl")
 
 function start_server(config_path)
     @info "Launching server... please wait"
 
-    warmup_cache(config_path)
-
     @info "Parsing configuration from $(config_path)..."
     config = TOML.parsefile(config_path)
 
+    @info "Initialising regional data and setting up tile cache"
+    initialise_data_with_cache(config)
+
     @info "Setting up auth middleware and router."
     auth = get_auth_router(config)
 
-    @info "Setting up criteria routes..."
-    setup_criteria_routes(config, auth)
-
-    @info "Setting up region routes..."
-    setup_region_routes(config, auth)
-
-    @info "Setting up tile routes..."
-    setup_tile_routes(config, auth)
-
-    @info "Setting up job routes..."
-    setup_job_routes(config, auth)
-
-    @info "Setting up admin routes..."
-    setup_admin_routes(config)
+    @info "Setting up utility routes..."
+    setup_utility_routes(config, auth)
 
     # Which host should we listen on?
     host = config["server_config"]["HOST"]
+
     # Which port should we listen on?
     port = 8000
 
@@ -110,35 +67,15 @@ This is a blocking operation until the worker times out.
 function start_worker()
     @info "Initializing worker from environment variables..."
     worker = create_worker_from_env()
-    @info "Starting worker loop from ReefGuideAPI.jl"
+    @info "Parsing TOML config"
+    config = TOML.parsefile(worker.config.config_path)
+    @info "Loading regional data"
+    initialise_data_with_cache(config)
+    @info "Starting worker loop from ReefGuideAPI.jl with $(Threads.nthreads()) threads."
     start(worker)
     @info "Worker closed itself..."
 end
 
-export
-    RegionalCriteria,
-    criteria_data_map
-
-# Methods to assess/identify deployment "plots" of reef.
-export
-    assess_reef_site,
-    identify_edge_aligned_sites,
-    filter_sites,
-    output_geojson
-
-# Geometry handling
-export
-    create_poly,
-    create_bbox,
-    port_buffer_mask,
-    meters_to_degrees,
-    polygon_to_lines
-
-# Raster->Index interactions (defunct?)
-export
-    valid_slope_lon_inds,
-    valid_slope_lat_inds,
-    valid_flat_lon_inds,
-    valid_flat_lat_inds
+export start_worker, start_server
 
 end

src/assessment_methods/apply_criteria.jl

@@ -0,0 +1,105 @@
+"""Methods to filter criteria bounds over rasters and lookup tables"""
+
+"""
+CriteriaBounds combine lookup information for a given criteria, bounds, and a
+rule (function) which enforces it for a given value
+"""
+struct CriteriaBounds{F<:Function}
+    "The field ID of the criteria"
+    name::Symbol
+    "min"
+    lower_bound::Float32
+    "max"
+    upper_bound::Float32
+    "A function which takes a value and returns if matches the criteria"
+    rule::F
+
+    function CriteriaBounds(name::S, lb::S, ub::S)::CriteriaBounds where {S<:String}
+        lower_bound::Float32 = parse(Float32, lb)
+        upper_bound::Float32 = parse(Float32, ub)
+        func = (x) -> lower_bound .<= x .<= upper_bound
+
+        return new{Function}(Symbol(name), lower_bound, upper_bound, func)
+    end
+
+    function CriteriaBounds(
+        name::String, lb::Float32, ub::Float32
+    )::CriteriaBounds
+        func = (x) -> lb .<= x .<= ub
+        return new{Function}(Symbol(name), lb, ub, func)
+    end
+end
+
+"""
+Apply thresholds for each criteria.
+
+# Arguments
+- `criteria_stack` : RasterStack of criteria data for a given region
+- `lookup` : Lookup dataframe for the region
+- `criteria_bounds` : A vector of CriteriaBounds which contains named criteria
+  with min/max ranges and a function to apply.
+
+# Returns
+BitMatrix indicating locations within desired thresholds
+"""
+function filter_raster_by_criteria(
+    criteria_stack::RasterStack,
+    lookup::DataFrame,
+    criteria_bounds::Vector{CriteriaBounds}
+)::Raster
+    # Result store
+    data = falses(size(criteria_stack))
+
+    # Apply criteria
+    res_lookup = trues(nrow(lookup))
+    for filter::CriteriaBounds in criteria_bounds
+        res_lookup .= res_lookup .& filter.rule(lookup[!, filter.name])
+    end
+
+    tmp = lookup[res_lookup, [:lon_idx, :lat_idx]]
+    data[CartesianIndex.(tmp.lon_idx, tmp.lat_idx)] .= true
+
+    res = Raster(criteria_stack.Depth; data=sparse(data), missingval=0)
+    return res
+end
+
+"""
+Filters the slope table (which contains raster param values too) by building a
+bit mask AND'd for all thresholds
+"""
+function filter_lookup_table_by_criteria(
+    slope_table::DataFrame,
+    ruleset::Vector{CriteriaBounds}
+)::DataFrame
+    slope_table.all_crit .= 1
+
+    for threshold in ruleset
+        slope_table.all_crit =
+            slope_table.all_crit .& threshold.rule(slope_table[!, threshold.name])
+    end
+
+    return slope_table[BitVector(slope_table.all_crit), :]
+end
+
+"""
+    lookup_df_from_raster(raster::Raster, threshold::Union{Int64,Float64})::DataFrame
+
+Build a look up table identifying all pixels in a raster that meet a suitability threshold.
+
+# Arguments
+- `raster` : Raster of regional data
+- `threshold` : Suitability threshold value (greater or equal than)
+
+# Returns
+DataFrame containing indices, lon and lat for each pixel that is intended for further
+analysis.
+"""
+function lookup_df_from_raster(raster::Raster, threshold::Union{Int64,Float64})::DataFrame
+    criteria_matches::SparseMatrixCSC{Bool,Int64} = sparse(falses(size(raster)))
+    Rasters.read!(raster .>= threshold, criteria_matches)
+    indices::Vector{CartesianIndex{2}} = findall(criteria_matches)
+    indices_lon::Vector{Float64} = lookup(raster, X)[first.(Tuple.(indices))]
+    indices_lat::Vector{Float64} = lookup(raster, Y)[last.(Tuple.(indices))]
+
+    return DataFrame(; indices=indices, lons=indices_lon, lats=indices_lat)
+end

src/assessment_methods/assessment_methods.jl

@@ -0,0 +1,21 @@
+using
+    Statistics,
+    Proj,
+    LibGEOS,
+    GeometryBasics,
+    CoordinateTransformations,
+    Rasters,
+    StaticArrays,
+    NearestNeighbors
+
+import ArchGDAL as AG
+import GeoInterface as GI
+import GeoInterface.Wrappers as GIWrap
+import GeometryOps as GO
+import GeoDataFrames as GDF
+
+include("apply_criteria.jl")
+include("best_fit_polygons.jl")
+include("common_functions.jl")
+include("geom_ops.jl")
+include("site_identification.jl")

src/site_assessment/best_fit_polygons.jl → src/assessment_methods/best_fit_polygons.jl

@@ -1,6 +1,5 @@
 """Geometry-based assessment methods."""
 
-using NearestNeighbors
 
 # Tabular data assessment methods
 
@@ -85,6 +84,7 @@ function assess_reef_site(
     best_poly[argmax(score)],
     maximum(qc_flag)
 end
+
 function assess_reef_site(
     rel_pix::DataFrame,
     rotated::Vector{GI.Wrappers.Polygon},
@@ -512,68 +512,6 @@ function find_optimal_site_alignment(
     )
 end
 
-# Raster based assessment methods
-
-# """
-#     assess_reef_site(
-#         rst::Union{Raster,RasterStack},
-#         geom::GI.Wrappers.Polygon,
-#         ruleset::Dict{Symbol,Function};
-#         degree_step::Float64=15.0,
-#         start_rot::Float64=0.0,
-#         n_per_side::Int64=1
-#     )::Tuple{Float64,Int64,GI.Wrappers.Polygon}
-
-# Assess given reef site for it's suitability score at different specified rotations around the
-# initial reef-edge rotation.
-
-# # Arguments
-# - `rst` : Raster of suitability scores.
-# - `geom` : Initial site polygon with no rotation applied.
-# - `ruleset` : Criteria ruleset to apply to `rst` pixels when assessing which pixels are suitable.
-# - `degree_step` : Degree value to vary each rotation by. Default = 15 degrees.
-# - `start_rot` : Initial rotation used to align the site polygon with the nearest reef edge. Default = 0 degrees.
-# - `n_per_side` : Number of times to rotate polygon on each side (clockwise and anticlockwise). Default = 2 rotations on each side.
-
-# # Returns
-# - Highest score identified with rotating polygons.
-# - The index of the highest scoring rotation.
-# - The polygon with the highest score out of the assessed rotated polygons.
-# """
-# function assess_reef_site(
-#     rst::Union{Raster,RasterStack},
-#     geom::GI.Wrappers.Polygon;
-#     degree_step::Float64=15.0,
-#     start_rot::Float64=0.0,
-#     n_per_side::Int64=2
-# )::Tuple{Float64,Int64,GI.Wrappers.Polygon}
-#     rotations =
-#         (start_rot - (degree_step * n_per_side)):degree_step:(start_rot + (degree_step * n_per_side))
-#     n_rotations = length(rotations)
-#     score = zeros(n_rotations)
-#     best_poly = Vector(undef, n_rotations)
-
-#     target_crs = convert(EPSG, GI.crs(rst))
-#     for (j, r) in enumerate(rotations)
-#         rot_geom = rotate_geom(geom, r, target_crs)
-#         c_rst = crop(rst; to=rot_geom)
-#         if !all(size(c_rst) .> (0, 0))
-#             @warn "No data found!"
-#             continue
-#         end
-
-#         score[j] = mean(c_rst)
-#         best_poly[j] = rot_geom
-
-#         if score[j] > 95.0
-#             # Found a good enough rotation
-#             break
-#         end
-#     end
-
-#     return score[argmax(score)], argmax(score) - (n_per_side + 1), best_poly[argmax(score)]
-# end
-
 """
     assess_reef_site(
         rst::Union{Raster,RasterStack},