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/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. diff --git a/src/criteria_assessment/tiles.jl b/src/criteria_assessment/tiles.jl index 1565e8d..6c275b4 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 = [ @@ -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,43 @@ 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 @@ -241,12 +212,17 @@ 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))