From c1e94c59b1ba110aeed287439b15a0134496d867 Mon Sep 17 00:00:00 2001 From: Takuya Iwanaga Date: Wed, 2 Oct 2024 18:11:00 +1000 Subject: [PATCH 1/5] Add Interpolations dependency --- Manifest.toml | 2 +- Project.toml | 1 + src/criteria_assessment/tiles.jl | 2 +- 3 files changed, 3 insertions(+), 2 deletions(-) diff --git a/Manifest.toml b/Manifest.toml index 5095d4b..0b173e6 100644 --- a/Manifest.toml +++ b/Manifest.toml @@ -2,7 +2,7 @@ julia_version = "1.10.5" manifest_format = "2.0" -project_hash = "6be7d23c68032d216194a57787647d02f3c8486f" +project_hash = "0416a1a4325488710cefe292c19281f43ae678ce" [[deps.AbstractFFTs]] deps = ["LinearAlgebra"] diff --git a/Project.toml b/Project.toml index 4ff203f..1eb9a0e 100644 --- a/Project.toml +++ b/Project.toml @@ -23,6 +23,7 @@ Glob = "c27321d9-0574-5035-807b-f59d2c89b15c" HTTP = "cd3eb016-35fb-5094-929b-558a96fad6f3" ImageIO = "82e4d734-157c-48bb-816b-45c225c6df19" Images = "916415d5-f1e6-5110-898d-aaa5f9f070e0" +Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" JSONWebTokens = "9b8beb19-0777-58c6-920b-28f749fee4d3" LibGEOS = "a90b1aa1-3769-5649-ba7e-abc5a9d163eb" diff --git a/src/criteria_assessment/tiles.jl b/src/criteria_assessment/tiles.jl index 1565e8d..1156734 100644 --- a/src/criteria_assessment/tiles.jl +++ b/src/criteria_assessment/tiles.jl @@ -2,7 +2,7 @@ Helper methods to support tiling """ -using ImageIO, Images +using ImageIO, Images, Interpolations # HTTP response headers for tile images const TILE_HEADERS = [ From ac83df92cf330fb1fe41455cc690d4c5146aa84e Mon Sep 17 00:00:00 2001 From: Takuya Iwanaga Date: Wed, 2 Oct 2024 18:11:57 +1000 Subject: [PATCH 2/5] Improved resampling, accounting for curvature --- src/criteria_assessment/tiles.jl | 94 +++++++++----------------------- 1 file changed, 25 insertions(+), 69 deletions(-) diff --git a/src/criteria_assessment/tiles.jl b/src/criteria_assessment/tiles.jl index 1156734..1502139 100644 --- a/src/criteria_assessment/tiles.jl +++ b/src/criteria_assessment/tiles.jl @@ -64,74 +64,29 @@ function _lon_lat_to_tile(zoom, lon, lat) return x, y end -""" - nearest(rst::Raster, tile_size::Tuple{Int, Int})::Matrix - -Resample a raster to a tile size using nearest neighbor interpolation. -This approach prioritising performance over accuracy. - -# Arguments -- `rst`: The input raster to be resampled. -- `tile_size`: The desired dimensions of the tile (lat, long). - -# Returns -Matrix with the resampled data. - -# Examples -```julia -large_raster = Raster(rand(UInt8, 14756, 14838); dims=(X(1:1:14756), Y(1:1:14838))) -small_matrix = nearest(large_raster, (256, 256)) -``` -""" -function nearest(rst::Raster, tile_size::Tuple{Int,Int})::Matrix - old_size = size(rst) - - # Important: must flip axes! - # Rasters.jl stores longitude along first dimension (rows) by default. - x_ratio = old_size[1] / tile_size[2] - y_ratio = old_size[2] / tile_size[1] - - resampled = zeros(eltype(rst), tile_size) - - Threads.@threads for lat in 1:tile_size[1] - for lon in 1:tile_size[2] - # Use area averaging for downsampling - x_start = max(1, floor(Int, (lon - 1) * x_ratio) + 1) - x_end = min(old_size[1], ceil(Int, lon * x_ratio)) - y_start = max(1, floor(Int, (lat - 1) * y_ratio) + 1) - y_end = min(old_size[2], ceil(Int, lat * y_ratio)) - - count_val = count(rst[x_start:x_end, y_start:y_end].data .> 0) - if count_val == 0 - continue - end - - sum_val = sum(rst[x_start:x_end, y_start:y_end].data) - resampled[lat, lon] = ceil(sum_val / count_val) - end - end - - return resampled -end +# Helper functions for Web Mercator projection +_lat_to_y(lat::Float64) = log(tan(π/4 + lat*π/360)) +_y_to_lat(y::Float64) = 360/π * atan(exp(y)) - 90 """ - masked_nearest(rst::Raster, z::Int, x::Int, y::Int, tile_size::Tuple{Int,Int}, orig_rst_size::Tuple{Int,Int})::Matrix + adjusted_nearest(rst::Raster, z::Int, x::Int, y::Int, tile_size::Tuple{Int,Int}, orig_rst_size::Tuple{Int,Int})::Matrix Resample a raster using nearest neighbor interpolation when the tile includes area outside where data exists (e.g., viewing the globe where the data may appear in a small corner of -the tile). This approach prioritising performance over accuracy. +the tile). This approach attempts to account for planetary curvature while still maintaining +some performance. # Arguments - `rst`: The input raster to be resampled. - `z`: Tile zoom level requested. - `x`: x coordinate for requested tile. - `y`: y coordinate for the requested tile. -- `tile_size`: The desired dimensions of the tile (lat, long). Note reversed order of coordinates. +- `tile_size`: The desired dimensions of the tile (long, lat). # Returns Matrix with the resampled data. """ -function masked_nearest( +function adjusted_nearest( rst::Raster, z::Int, x::Int, @@ -144,27 +99,28 @@ function masked_nearest( # Bounds for the area of interest (AOI; where we have data) ((aoi_lon_min, aoi_lon_max), (aoi_lat_min, aoi_lat_max)) = Rasters.bounds(rst) - # Create an empty tile (lat/long) - lat_size = tile_size[1] - long_size = tile_size[2] - tile = fill(0.0, lat_size, long_size) + # Create an empty tile (long/lat) + long_size, lat_size = tile_size + tile = zeros(lat_size, long_size) + + # Generate longitude and latitude arrays for the tile + lons = @. mod(t_lon_min + (t_lon_max - t_lon_min) * ((1:long_size) - 1) / (long_size - 1) + 180, 360) - 180 + lats = @. _y_to_lat(_lat_to_y(t_lat_max) - (_lat_to_y(t_lat_max) - _lat_to_y(t_lat_min)) * ((1:lat_size) - 1) / (lat_size - 1)) - lons = @. t_lon_min + (t_lon_max - t_lon_min) * ((1:long_size) - 1) / (long_size - 1) - lats = @. t_lat_max - (t_lat_max - t_lat_min) * ((1:lat_size) - 1) / (lat_size - 1) + # Determine which points are within the area of interest in_lons = aoi_lon_min .<= lons .<= aoi_lon_max in_lats = aoi_lat_min .<= lats .<= aoi_lat_max # Sample data that is within area of interest - data_x = - round.( - Int, (lons[in_lons] .- aoi_lon_min) / (aoi_lon_max - aoi_lon_min) * size(rst, 1) - ) - data_y = - round.( - Int, (aoi_lat_max .- lats[in_lats]) / (aoi_lat_max - aoi_lat_min) * size(rst, 2) - ) - - tile[in_lats, in_lons] .= rst[data_x, data_y]' + for (i, lon) in enumerate(lons), (j, lat) in enumerate(lats) + if in_lons[i] && in_lats[j] + x_idx = round(Int, (lon - aoi_lon_min) / (aoi_lon_max - aoi_lon_min) * (size(rst, 1) - 1)) + 1 + y_idx = round(Int, (_lat_to_y(aoi_lat_max) - _lat_to_y(lat)) / (_lat_to_y(aoi_lat_max) - _lat_to_y(aoi_lat_min)) * (size(rst, 2) - 1)) + 1 + x_idx = clamp(x_idx, 1, size(rst, 1)) + y_idx = clamp(y_idx, 1, size(rst, 2)) + tile[j, i] = rst[x_idx, y_idx] + end + end return tile end From d7e43ced6dd4569157edb02aa4fa0bcc6cdee90c Mon Sep 17 00:00:00 2001 From: Takuya Iwanaga Date: Wed, 2 Oct 2024 18:12:31 +1000 Subject: [PATCH 3/5] Signpost why additional params are needed --- src/criteria_assessment/query_thresholds.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/criteria_assessment/query_thresholds.jl b/src/criteria_assessment/query_thresholds.jl index dd8bb95..2fffaf2 100644 --- a/src/criteria_assessment/query_thresholds.jl +++ b/src/criteria_assessment/query_thresholds.jl @@ -206,8 +206,8 @@ applied to a set of criteria. - `reg_criteria` : RegionalCriteria to assess - `rtype` : reef type to assess (`:slopes` or `:flats`) - `crit_map` : List of criteria thresholds to apply (see `apply_criteria_thresholds()`) -- `lons` : Longitudinal extent (min and max) -- `lats` : Latitudinal extent (min and max) +- `lons` : Longitudinal extent (min and max, required when generating masks for tiles) +- `lats` : Latitudinal extent (min and max, required when generating masks for tiles) # Returns True/false mask indicating locations within desired thresholds. From 69be1593fd98e23ac8df6a9a076eae54e7ce6cb7 Mon Sep 17 00:00:00 2001 From: Takuya Iwanaga Date: Wed, 2 Oct 2024 18:13:44 +1000 Subject: [PATCH 4/5] Adjust when direct resampling is applied The further user is zoomed in, the less need there is to perform additional calculations to account for the curvature of the Earth, so can do a direct resampling. --- src/criteria_assessment/tiles.jl | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/criteria_assessment/tiles.jl b/src/criteria_assessment/tiles.jl index 1502139..f0dbd34 100644 --- a/src/criteria_assessment/tiles.jl +++ b/src/criteria_assessment/tiles.jl @@ -197,12 +197,15 @@ function setup_tile_routes(config, auth) # Using if block to avoid type instability @debug "Thread $(thread_id) - $(now()) : Creating PNG (with transparency)" if any(size(mask_data) .> tile_size(config)) - if any(size(mask_data) .== size(reg_assess_data[reg].stack)) || (z < 8) + if any(size(mask_data) .== size(reg_assess_data[reg].stack)) || (z < 12) # Account for geographic positioning when zoomed out further than # raster area - resampled = masked_nearest(mask_data, z, x, y, tile_size(config)) + resampled = adjusted_nearest(mask_data, z, x, y, tile_size(config)) else - resampled = nearest(mask_data, tile_size(config)) + # Zoomed in close so less need to account for curvature + # BSpline(Constant()) is equivalent to nearest neighbor. + # See details in: https://juliaimages.org/ImageTransformations.jl/stable/reference/#Low-level-warping-API + resampled = imresize(mask_data, tile_size(config); method=BSpline(Constant())) end img = zeros(RGBA, size(resampled)) From 85b4d7889694acb146ae3033fd7a0331ebaa93ce Mon Sep 17 00:00:00 2001 From: Takuya Iwanaga Date: Wed, 2 Oct 2024 18:30:00 +1000 Subject: [PATCH 5/5] Appease the formatter --- src/criteria_assessment/tiles.jl | 31 ++++++++++++++++++++++++------- 1 file changed, 24 insertions(+), 7 deletions(-) diff --git a/src/criteria_assessment/tiles.jl b/src/criteria_assessment/tiles.jl index f0dbd34..6c275b4 100644 --- a/src/criteria_assessment/tiles.jl +++ b/src/criteria_assessment/tiles.jl @@ -65,8 +65,8 @@ function _lon_lat_to_tile(zoom, lon, lat) end # Helper functions for Web Mercator projection -_lat_to_y(lat::Float64) = log(tan(π/4 + lat*π/360)) -_y_to_lat(y::Float64) = 360/π * atan(exp(y)) - 90 +_lat_to_y(lat::Float64) = log(tan(π / 4 + lat * π / 360)) +_y_to_lat(y::Float64) = 360 / π * atan(exp(y)) - 90 """ adjusted_nearest(rst::Raster, z::Int, x::Int, y::Int, tile_size::Tuple{Int,Int}, orig_rst_size::Tuple{Int,Int})::Matrix @@ -104,8 +104,14 @@ function adjusted_nearest( tile = zeros(lat_size, long_size) # Generate longitude and latitude arrays for the tile - lons = @. mod(t_lon_min + (t_lon_max - t_lon_min) * ((1:long_size) - 1) / (long_size - 1) + 180, 360) - 180 - lats = @. _y_to_lat(_lat_to_y(t_lat_max) - (_lat_to_y(t_lat_max) - _lat_to_y(t_lat_min)) * ((1:lat_size) - 1) / (lat_size - 1)) + lons = @. mod( + t_lon_min + (t_lon_max - t_lon_min) * ((1:long_size) - 1) / (long_size - 1) + 180, + 360 + ) - 180 + lats = @. _y_to_lat( + _lat_to_y(t_lat_max) - + (_lat_to_y(t_lat_max) - _lat_to_y(t_lat_min)) * ((1:lat_size) - 1) / (lat_size - 1) + ) # Determine which points are within the area of interest in_lons = aoi_lon_min .<= lons .<= aoi_lon_max @@ -114,8 +120,17 @@ function adjusted_nearest( # Sample data that is within area of interest for (i, lon) in enumerate(lons), (j, lat) in enumerate(lats) if in_lons[i] && in_lats[j] - x_idx = round(Int, (lon - aoi_lon_min) / (aoi_lon_max - aoi_lon_min) * (size(rst, 1) - 1)) + 1 - y_idx = round(Int, (_lat_to_y(aoi_lat_max) - _lat_to_y(lat)) / (_lat_to_y(aoi_lat_max) - _lat_to_y(aoi_lat_min)) * (size(rst, 2) - 1)) + 1 + x_idx = + round( + Int, + (lon - aoi_lon_min) / (aoi_lon_max - aoi_lon_min) * (size(rst, 1) - 1) + ) + 1 + y_idx = + round( + Int, + (_lat_to_y(aoi_lat_max) - _lat_to_y(lat)) / + (_lat_to_y(aoi_lat_max) - _lat_to_y(aoi_lat_min)) * (size(rst, 2) - 1) + ) + 1 x_idx = clamp(x_idx, 1, size(rst, 1)) y_idx = clamp(y_idx, 1, size(rst, 2)) tile[j, i] = rst[x_idx, y_idx] @@ -205,7 +220,9 @@ function setup_tile_routes(config, auth) # Zoomed in close so less need to account for curvature # BSpline(Constant()) is equivalent to nearest neighbor. # See details in: https://juliaimages.org/ImageTransformations.jl/stable/reference/#Low-level-warping-API - resampled = imresize(mask_data, tile_size(config); method=BSpline(Constant())) + resampled = imresize( + mask_data, tile_size(config); method=BSpline(Constant()) + ) end img = zeros(RGBA, size(resampled))