Merge pull request #7 from open-AIMS/geojson_output

Add GeoJSON file outputs for site polygons.

src/ReefGuideAPI.jl

@@ -244,7 +244,8 @@ export
 export
     assess_reef_site,
     identify_potential_sites_edges,
-    filter_intersecting_sites
+    filter_sites,
+    output_geojson
 
 # Geometry handling
 export

src/geom_handlers/common_assessment.jl

@@ -1,10 +1,10 @@
-""" Functions common to both site_assessment methods. """
+""" Functions common to both site_assessment methods."""
 
+using Dates
 using LinearAlgebra
 
 include("geom_ops.jl")
 
-# Additional functions for reef-edge alignment processing.
 """
     meters_to_degrees(x, lat)
 
@@ -73,12 +73,18 @@ function line_angle(a::T, b::T)::Float64 where {T<:Vector{Tuple{Float64,Float64}
 end
 
 """
-    filter_far_polygons(gdf::DataFrame, pixel::GIWrap.Point, lat::Float64)::BitVector
+    filter_far_polygons(gdf::DataFrame, pixel::GIWrap.Point, lat::Float64, dist::Union{Int64,Float64})::BitVector
 
 Filter out reefs that are > 10km from the target pixel (currently hardcoded threshold).
 """
-function filter_far_polygons(gdf::DataFrame, pixel::GeometryBasics.Point, lat::Float64)::BitVector
-    return (GO.distance.(GO.centroid.(gdf[:, first(GI.geometrycolumns(gdf))]), [pixel]) .< meters_to_degrees(10000, lat))
+function filter_far_polygons(
+    gdf::DataFrame,
+    pixel::GeometryBasics.Point,
+    lat::Float64,
+    dist::Union{Int64,Float64}
+)::BitVector
+    geoms = gdf[:, first(GI.geometrycolumns(gdf))]
+    return (GO.distance.(GO.centroid.(geoms), [pixel]) .< meters_to_degrees(dist, lat))
 end
 
 """
@@ -101,7 +107,7 @@ and is buffered by `res` distance.
 - `res` : Buffer distance (resolution of input raster search data).
 
 # Returns
-- Initial search box geometry.
+Initial search box geometry.
 """
 function initial_search_box(
     (lon, lat),
@@ -122,7 +128,7 @@ function initial_search_box(
 end
 
 """
-    identify_closest_edge(
+    closest_reef_edge(
         pixel::GeometryBasics.Point{2, Float64},
         reef_lines::Vector{GeometryBasics.Line{2, Float64}}
     )::Vector{Tuple{Float64, Float64}}
@@ -132,8 +138,11 @@ Find the nearest line in `reef_lines` to a point `pixel`.
 # Arguments
 - `pixel` : Target point geometry.
 - `reef_lines` : Vector containing lines for comparison.
+
+# Returns
+Coordinates of the reef edge line that is closest to the target `pixel`. Returned in Tuples.
 """
-function identify_closest_edge(
+function closest_reef_edge(
     pixel::GeometryBasics.Point{2,Float64},
     reef_lines::Vector{GeometryBasics.Line{2,Float64}}
 )::Vector{Tuple{Float64,Float64}}
@@ -147,7 +156,8 @@ end
         pixel::GeometryBasics.Point{2, Float64},
         geom_buff::GI.Wrappers.Polygon,
         gdf::DataFrame,
-        reef_outlines::Vector{Vector{GeometryBasics.Line{2, Float64}}}
+        reef_outlines::Vector{Vector{GeometryBasics.Line{2, Float64}}};
+        search_buffer::Union{Int64,Float64}=20000.0
     )::Float64
 
 Identifies the closest edge to the target `pixel`/'`geom_buff` and returns the initial rotation
@@ -163,9 +173,10 @@ function initial_search_rotation(
     pixel::GeometryBasics.Point{2,Float64},
     geom_buff::GI.Wrappers.Polygon,
     gdf::DataFrame,
-    reef_outlines::Vector{Vector{GeometryBasics.Line{2,Float64}}}
+    reef_outlines::Vector{Vector{GeometryBasics.Line{2,Float64}}};
+    search_buffer::Union{Int64,Float64}=20000.0
 )::Float64
-    distance_indices = filter_far_polygons(gdf, pixel, pixel[2])
+    distance_indices = filter_far_polygons(gdf, pixel, pixel[2], search_buffer)
     reef_lines = reef_outlines[distance_indices]
     reef_lines = reef_lines[
         GO.within.([pixel], gdf[distance_indices, first(GI.geometrycolumns(gdf))])
@@ -183,7 +194,7 @@ function initial_search_rotation(
         reef_lines = vcat(reef_lines...)
     end
 
-    edge_line = identify_closest_edge(pixel, reef_lines)
+    edge_line = closest_reef_edge(pixel, reef_lines)
 
     # Calculate the angle between the two lines
     edge_bearing = line_angle([(0.0, 5.0), (0.0, 0.0)], from_zero(edge_line))
@@ -196,16 +207,21 @@ function initial_search_rotation(
 end
 
 """
-    filter_intersecting_sites(res_df::DataFrame)::DataFrame
+    filter_sites(res_df::DataFrame)::DataFrame
 
 Filter out sites where the qc_flag indicates a suitabiltiy < `surr_threshold` in searching.
 Identify and keep the highest scoring site polygon where site polygons are overlapping.
 
 # Arguments
-- `res_df` : Results DataFrame containing potential site polygons (output from
-`identify_potential_sites()` or `identify_potential_sites_edges()`).
+- `res_df` : Results DataFrame containing potential site polygons
+             (output from `identify_potential_sites()` or `identify_potential_sites_edges()`).
+
+# Returns
+DataFrame containing only the highest scoring sites where site polygons intersect, and
+containing only sites with scores greater than the `surr_threshold` specified in
+`identify_potential_sites_edges()` (default=0.33).
 """
-function filter_intersecting_sites(res_df::DataFrame)::DataFrame
+function filter_sites(res_df::DataFrame)::DataFrame
     res_df.row_ID = 1:size(res_df, 1)
     ignore_list = []
 
@@ -231,5 +247,65 @@ function filter_intersecting_sites(res_df::DataFrame)::DataFrame
         end
     end
 
-    return res_df[res_df.row_ID .∉ [unique(ignore_list)], :]
+    rename!(res_df, :poly => :geometry)
+
+    return res_df[res_df.row_ID .∉ [unique(ignore_list)], Not(:qc_flag, :row_ID)]
+end
+
+"""
+    output_geojson(
+        df::DataFrame,
+        region::String,
+        output_dir::String
+    )::Nothing
+
+Writes out GeoJSON file to a target directory. Output file will be located at location:
+"`output_dir`/output_sites_`region`.geojson"
+
+# Arguments
+- `df` : DataFrame intended for writing to geojson file.
+- `region` : Region name for labelling output file.
+- `output_dir` : Directory to write geojson file to.
+"""
+function output_geojson(
+    df::DataFrame,
+    region::String,
+    output_dir::String
+)::Nothing
+    GDF.write(
+        joinpath(
+            output_dir,
+            "output_sites_$(region).geojson"
+        ),
+        df;
+        crs=GI.crs(first(df.geometry))
+    )
+
+    return nothing
+end
+
+"""
+    identify_search_pixels(input_raster::Raster, criteria_function)::DataFrame
+
+Identifies all pixels in an input raster that return true for the function `criteria_function`.
+
+# Arguments
+- `input_raster` : Raster containing pixels for the target region.
+- `criteria_function` : Function that returns a boolean value for each pixel in `input_raster`. Pixels that return true will be targetted in analysis.
+
+# Returns
+DataFrame containing indices, lon and lat for each pixel that is intended for further analysis.
+"""
+function identify_search_pixels(input_raster::Raster, criteria_function)::DataFrame
+    pixels = trim(criteria_function(scan_locs))
+    indices = findall(pixels)
+    indices_lon = Vector{Union{Missing, Float64}}(missing, size(indices, 1))
+    indices_lat = Vector{Union{Missing, Float64}}(missing, size(indices, 1))
+
+    for (j, index) in enumerate(indices)
+        indices_lon[j] = dims(pixels, X)[index[1]]
+        indices_lat[j] = dims(pixels, Y)[index[2]]
+    end
+
+    return DataFrame(indices = indices, lon = indices_lon, lat = indices_lat)
 end

src/geom_handlers/geom_ops.jl

@@ -87,7 +87,11 @@ function get_points(geom)
     end
 end
 
-function rotate_geom(geom, degrees::Float64)
+function rotate_geom(
+    geom,
+    degrees::Float64,
+    target_crs::GeoFormatTypes.CoordinateReferenceSystemFormat
+)
     degrees == 0.0 && return geom
 
     theta = deg2rad(degrees)
@@ -113,7 +117,7 @@ function rotate_geom(geom, degrees::Float64)
         new_points[i] = rotate_point(new_points[i])
     end
 
-    return create_poly(new_points, GI.crs(geom))
+    return create_poly(new_points, target_crs)
 end
 
 """

src/geom_handlers/parq_site_assessment.jl

@@ -7,22 +7,38 @@ include("common_assessment.jl")
     assess_reef_site(
         rel_pix::DataFrame,
         geom::GI.Wrappers.Polygon,
-        max_count::Float64;
+        max_count::Float64,
+        target_crs::GeoFormatTypes.CoordinateReferenceSystemFormat;
         degree_step::Float64=15.0,
         start_rot::Float64=0.0,
         n_per_side::Int64=2,
         surr_threshold::Float64=0.33
     )::Tuple{Float64,Int64,GI.Wrappers.Polygon,Int64}
 
-Assesses the rotations of a search box `geom` for their pixel score. Returns the highest score,
-rotation step, polygon and a quality control flag for each assessment. Rotation steps are returned
-so that the `start_rot` angle is 0, rotations anti-clockwise are negative and rotations clockwise are
+Assesses the rotations of a search box `geom` for their suitability score (calculated as the
+proportion of pixels that meet all specified criteria thresholds). Search box rotation steps
+are returned so that the `start_rot` angle is 0, rotations anti-clockwise are negative and
+rotations clockwise are
 positive.
+
+# Arguments
+- `rel_pix` : DataFrame containing the point data for pixels that are within maxmimum user search box dimensions from a pixel.
+- `geom` : Starting search box for assessment.
+- `max_count` : The maximum number of pixels that can intersect the search box (used to standardise scores between 0 and 1).
+- `target_crs` : Coordinate Reference System used for analysis vector and raster data.
+- `degree_step` : Step to vary the search box rotations.
+- `start_rot` : Starting angle rotation that aligns the box with the closest reef edge.
+- `n_per_side` : Number of rotations to perform around the starting search box angle.
+- `surr_threshold` : Suitability threshold, below which sites are excluded from result sets.
+
+# Returns
+Returns the highest score, rotation step, polygon and a quality control flag for each assessment.
 """
 function assess_reef_site(
     rel_pix::DataFrame,
     geom::GI.Wrappers.Polygon,
-    max_count::Float64;
+    max_count::Float64,
+    target_crs::GeoFormatTypes.CoordinateReferenceSystemFormat;
     degree_step::Float64=15.0,
     start_rot::Float64=0.0,
     n_per_side::Int64=2,
@@ -35,7 +51,7 @@ function assess_reef_site(
     qc_flag = zeros(Int64, n_rotations)
 
     for (j, r) in enumerate(rotations)
-        rot_geom = rotate_geom(geom, r)
+        rot_geom = rotate_geom(geom, r, target_crs)
         score[j] = size(rel_pix[GO.intersects.([rot_geom], rel_pix.geometry), :], 1) / max_count
         best_poly[j] = rot_geom
 
@@ -51,8 +67,8 @@ end
 """
     identify_potential_sites_edges(
         df::DataFrame,
-        indices_pixels::Raster,
-        indices::Vector{CartesianIndex{2}},
+        search_pixels::DataFrame,
+        res::Float64,
         gdf::DataFrame,
         x_dist::Union{Int64,Float64},
         y_dist::Union{Int64,Float64},
@@ -64,16 +80,16 @@ end
         surr_threshold::Float64=0.33
     )::DataFrame
 
-Identify the most suitable site polygons for each pixel in the `search_pixels` raster where
-`indices` denotes which pixels to check for suitability. `x_dist` and `y_dist` are x and y
-lengths of the search polygon in meters. A buffer of `rst_stack` resolution is applied to the search box.
-And angle from a pixel to a reef edge is identified and used for searching with custom rotation
-parameters. Method is currently opperating for CRS in degrees units.
+Identify the most suitable site polygons for each pixel in the `search_pixels` DataFrame.
+`x_dist` and `y_dist` are x and y lengths of the search polygon in meters. A buffer of the
+raster files' resolution is applied to the search box. And angle from a pixel to a reef edge
+is identified and used for searching with custom rotation parameters.
+Method is currently opperating for CRS in degrees units.
 
 # Arguments
 - `df` : DataFrame containing environmental variables for assessment.
-- `indices_pixels` : Raster that matches indices for lon/lat information.
-- `indices` : Vector of CartesianIndices noting pixels to assess sites.
+- `search_pixels` : DataFrame containing lon and lat columns for each pixel that is intended for analysis.
+- `res` : Resolution of the original raster pixels. Can by found via `abs(step(dims(raster, X)))`.
 - `gdf` : GeoDataFrame containing the reef outlines used to align the search box edge.
 - `x_dist` : Length of horizontal side of search box (in meters).
 - `y_dist` : Length of vertical side of search box (in meters).
@@ -89,8 +105,8 @@ DataFrame containing highest score, rotation and polygon for each assessment at
 """
 function identify_potential_sites_edges(
     df::DataFrame,
-    indices_pixels::Raster,
-    indices::Vector{CartesianIndex{2}},
+    search_pixels::DataFrame,
+    res::Float64,
     gdf::DataFrame,
     x_dist::Union{Int64,Float64},
     y_dist::Union{Int64,Float64},
@@ -103,21 +119,20 @@ function identify_potential_sites_edges(
 )::DataFrame
     reef_lines = reef_lines[gdf.management_area .== region]
     gdf = gdf[gdf.management_area .== region, :]
-    res = abs(step(dims(indices_pixels, X)))
     max_count = (
         (x_dist / degrees_to_meters(res, mean(indices_pixels.dims[2]))) *
         ((y_dist + 2 * degrees_to_meters(res, mean(indices_pixels.dims[2]))) /
          degrees_to_meters(res, mean(indices_pixels.dims[2])))
     )
 
     # Search each location to assess
-    best_score = zeros(length(indices))
-    best_poly = Vector(undef, length(indices))
-    best_rotation = zeros(Int64, length(indices))
-    quality_flag = zeros(Int64, length(indices))
-    for (i, index) in enumerate(indices)
-        lon = dims(indices_pixels, X)[index[1]]
-        lat = dims(indices_pixels, Y)[index[2]]
+    best_score = zeros(length(search_pixels.lon))
+    best_poly = Vector(undef, length(search_pixels.lon))
+    best_rotation = zeros(Int64, length(search_pixels.lon))
+    quality_flag = zeros(Int64, length(search_pixels.lon))
+    for (i, index) in enumerate(eachrow(search_pixels))
+        lon = index.lon
+        lat = index.lat
         geom_buff = initial_search_box((lon, lat), x_dist, y_dist, target_crs, res)
 
         pixel = GO.Point(lon, lat)
@@ -136,7 +151,8 @@ function identify_potential_sites_edges(
         b_score, b_rot, b_poly, qc_flag = assess_reef_site(
             rel_pix,
             geom_buff,
-            max_count;
+            max_count,
+            target_crs;
             degree_step=degree_step,
             start_rot=rot_angle,
             n_per_side=n_rot_p_side,