From c240c89ba74107a8afef430c8434508d40747252 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Timoth=C3=A9e=20Poisot?= <timothee.poisot@umontreal.ca>
Date: Mon, 28 Jul 2025 12:02:10 -0400
Subject: [PATCH 01/19] feat(layers): quantile transfer

---
 SimpleSDMLayers/src/SimpleSDMLayers.jl  |  7 +++++++
 SimpleSDMLayers/src/quantiletransfer.jl | 24 ++++++++++++++++++++++++
 2 files changed, 31 insertions(+)
 create mode 100644 SimpleSDMLayers/src/quantiletransfer.jl

diff --git a/SimpleSDMLayers/src/SimpleSDMLayers.jl b/SimpleSDMLayers/src/SimpleSDMLayers.jl
index de9647e5e..8c068481f 100644
--- a/SimpleSDMLayers/src/SimpleSDMLayers.jl
+++ b/SimpleSDMLayers/src/SimpleSDMLayers.jl
@@ -9,6 +9,7 @@ rasters.
 """
 module SimpleSDMLayers
 
+using LinearAlgebra
 using TestItems
 
 import ArchGDAL
@@ -86,4 +87,10 @@ export cellarea
 include("reclassify.jl")
 export reclassify
 
+# Quantile transfer
+include("quantiletransfer.jl")
+export quantiletransfer, quantiletransfer!
+
+# Shift and rotate
+
 end # module
diff --git a/SimpleSDMLayers/src/quantiletransfer.jl b/SimpleSDMLayers/src/quantiletransfer.jl
new file mode 100644
index 000000000..7c27777ed
--- /dev/null
+++ b/SimpleSDMLayers/src/quantiletransfer.jl
@@ -0,0 +1,24 @@
+"""
+    quantiletransfer!(target::SDMLayer{T}, reference::SDMLayer; n::Int=10) where {T <: Real}
+
+Replaces the values in the `target` layer so that they follow the distribution
+of the values in the `reference` layer. This works by (i) identifying the
+quantile of each cell in the target layer, then (ii) replacing this with the
+value for the same quantile in the reference layer.
+"""
+function quantiletransfer!(target::SDMLayer{T}, reference::SDMLayer; n::Int=10) where {T <: Real}
+    Qt = quantize(target, round(Int, n * sqrt(count(reference))))
+    target.grid = quantile(reference, Qt.grid)
+    return target
+end
+
+"""
+    quantiletransfer(target, reference)
+
+Non-mutating version of `quantiletransfer!`
+"""
+function quantiletransfer(target, reference)
+    t = deepcopy(target)
+    quantiletransfer!(t, reference)
+    return t
+end
\ No newline at end of file

From 949516d1ff8a75994d1fe9157b788e72c789c25e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Timoth=C3=A9e=20Poisot?= <timothee.poisot@umontreal.ca>
Date: Mon, 28 Jul 2025 12:02:18 -0400
Subject: [PATCH 02/19] semver(layers): v1.4.0

---
 SimpleSDMLayers/CHANGELOG.md | 5 +++++
 SimpleSDMLayers/Project.toml | 4 +++-
 2 files changed, 8 insertions(+), 1 deletion(-)

diff --git a/SimpleSDMLayers/CHANGELOG.md b/SimpleSDMLayers/CHANGELOG.md
index 3c773a6ac..dc6cf0794 100644
--- a/SimpleSDMLayers/CHANGELOG.md
+++ b/SimpleSDMLayers/CHANGELOG.md
@@ -5,6 +5,11 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
+## `v1.4.0`
+
+- **added** functions to support shift and rotate
+- **added** a function for quantile transfer
+
 ## `v1.3.1`
 
 - **fixed** plotting support to use `Makie` instead of `MakieCore`
diff --git a/SimpleSDMLayers/Project.toml b/SimpleSDMLayers/Project.toml
index 7eeefe4f8..b84d02d04 100644
--- a/SimpleSDMLayers/Project.toml
+++ b/SimpleSDMLayers/Project.toml
@@ -1,11 +1,12 @@
 name = "SimpleSDMLayers"
 uuid = "2c645270-77db-11e9-22c3-0f302a89c64c"
 authors = ["Timothée Poisot <timothee.poisot@umontreal.ca>", "Gabriel Dansereau <gabriel.dansereau@umontreal.ca>"]
-version = "1.3.1"
+version = "1.4.0"
 
 [deps]
 ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"
 Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Proj = "c94c279d-25a6-4763-9509-64d165bea63e"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
@@ -28,6 +29,7 @@ NeutralLandscapesExtension = "NeutralLandscapes"
 ArchGDAL = "0.10"
 Clustering = "0.15"
 Distances = "0.10"
+LinearAlgebra = "1.11.0"
 Makie = "0.21, 0.22, 0.23"
 MultivariateStats = "0.10"
 NeutralLandscapes = "0.1"

From 4ae5b6928c77ae5f554d43faf1c0c5e302e2ec41 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Timoth=C3=A9e=20Poisot?= <timothee.poisot@umontreal.ca>
Date: Mon, 28 Jul 2025 12:02:27 -0400
Subject: [PATCH 03/19] feat(layers)?: shift and rotate

---
 SimpleSDMLayers/src/shift-and-rotate.jl | 102 ++++++++++++++++++++++++
 1 file changed, 102 insertions(+)
 create mode 100644 SimpleSDMLayers/src/shift-and-rotate.jl

diff --git a/SimpleSDMLayers/src/shift-and-rotate.jl b/SimpleSDMLayers/src/shift-and-rotate.jl
new file mode 100644
index 000000000..50537868b
--- /dev/null
+++ b/SimpleSDMLayers/src/shift-and-rotate.jl
@@ -0,0 +1,102 @@
+"""
+    _rotate_coordinate(lon, lat, roll, pitch, yaw)
+
+Transforms a coordinate under a rotation of the Earth, represented as a unit
+sphere. The rotation is specified using the roll (α, x-axis rotation), pitch (β,
+y-axis rotation), and yaw (γ, z-axis rotation), all angles are given in degrees.
+The returned point is given as a pair of lon, lat coordinates.
+"""
+function _rotate_coordinate(lon, lat, roll, pitch, yaw)
+
+    # Move everything to rads
+    α = deg2rad(roll)
+    β = deg2rad(pitch)
+    γ = deg2rad(yaw)
+
+    # Coordinates to unit sphere position
+    x = cos(deg2rad(lat)) * cos(deg2rad(lon))
+    y = cost(deg2rad(lat)) * sin(deg2rad(lon))
+    z = sin(deg2rad(lat))
+
+    # Point vector
+    p = [x, y, z]
+
+    # Rotation matrix for the roll
+    Rx = [1 0 0; 0 cos(α) -sin(α); 0 sin(α) cos(α)]
+    # Rotation matrix for the pitch
+    Ry = [cos(β) 0 sin(β); 0 1 0; -sin(β) - cos(β)]
+    # Rotation matrix for the yaw
+    Rz = [cos(γ) -sin(γ) 0; sin(γ) cos(γ) 0; 0 0 1]
+
+    # Rotation matrix (x, then y, then z)
+    R = Rz * Ry * Rx
+
+    # New point
+    X, Y, Z = R * p
+
+    # New longitude and latitude
+    projected = (atan(Y, X), asin(Z))
+    
+    # Return the new point in degrees
+    return rad2deg.(projected)
+end
+
+#=
+
+function generate_valid_rotation(reference, pool; roll=(-5.0, 5.0), pitch=(-5.0, 5.0), yaw=(-5.0, 5.0))
+    E, N, K = eastings(reference), northings(reference), keys(reference)
+    still_searching = true
+    Θ = (0.0, 0.0, 0.0)
+    while still_searching
+        α = rand() * (roll[2]-roll[1]) + roll[1] 
+        β = rand() * (pitch[2]-pitch[1]) + pitch[1]
+        γ = rand() * (yaw[2]-yaw[1]) + yaw[1]
+        Θ = (α, β, γ)
+        still_searching = false
+        for k in K
+            lat, lon = rotate_earth_coordinates(N[k[1]], E[k[2]], Θ...)
+            try
+                nv = pool[lon, lat]
+                if isnothing(nv)
+                    still_searching = true
+                    break
+                end
+            catch err
+                still_searching = true
+                break
+            end
+        end
+    end
+    return Θ
+end
+=#
+
+#=
+
+temp = SDMLayer(RasterData(WorldClim2, BioClim); layer=1, resolution=2.5)
+pol = getpolygon(PolygonData(NaturalEarth, Countries))["Czech Republic"]
+ref = trim(mask(temp, pol))
+lat, lon = northings(ref), eastings(ref)
+
+# rotation, up/down, left/right
+angles = generate_valid_rotation(ref, temp)
+@info angles
+
+sar = deepcopy(ref)
+
+for k in keys(ref)
+    coords = (lat[k[1]], lon[k[2]])
+    nlat, nlon = rotate_earth_coordinates(coords..., angles...)
+    sar.grid[k] = temp[nlon, nlat]
+end
+
+replaced = quantiletransfer(sar, ref)
+
+f = Figure()
+ax = Axis(f[1,1], aspect=DataAspect())
+ax2 = Axis(f[2,1])
+hm = heatmap!(ax, replaced, colormap=:navia)
+hist!(ax2, replaced, bins=100)
+Colorbar(f[1:2,2], hm)
+f
+=#
\ No newline at end of file

From 843516ad62b502b1bb9735313b9fad3ba59ebc0d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Timoth=C3=A9e=20Poisot?= <timothee.poisot@umontreal.ca>
Date: Mon, 28 Jul 2025 14:20:51 -0400
Subject: [PATCH 04/19] feat(layers)?: shift and rotate

---
 SimpleSDMLayers/src/SimpleSDMLayers.jl  |   2 +
 SimpleSDMLayers/src/shift-and-rotate.jl | 169 +++++++++++-------------
 2 files changed, 79 insertions(+), 92 deletions(-)

diff --git a/SimpleSDMLayers/src/SimpleSDMLayers.jl b/SimpleSDMLayers/src/SimpleSDMLayers.jl
index 8c068481f..1bc8ef6ef 100644
--- a/SimpleSDMLayers/src/SimpleSDMLayers.jl
+++ b/SimpleSDMLayers/src/SimpleSDMLayers.jl
@@ -10,6 +10,7 @@ rasters.
 module SimpleSDMLayers
 
 using LinearAlgebra
+using Statistics
 using TestItems
 
 import ArchGDAL
@@ -92,5 +93,6 @@ include("quantiletransfer.jl")
 export quantiletransfer, quantiletransfer!
 
 # Shift and rotate
+include("shift-and-rotate.jl")
 
 end # module
diff --git a/SimpleSDMLayers/src/shift-and-rotate.jl b/SimpleSDMLayers/src/shift-and-rotate.jl
index 50537868b..c05be48ae 100644
--- a/SimpleSDMLayers/src/shift-and-rotate.jl
+++ b/SimpleSDMLayers/src/shift-and-rotate.jl
@@ -1,102 +1,87 @@
-"""
-    _rotate_coordinate(lon, lat, roll, pitch, yaw)
-
-Transforms a coordinate under a rotation of the Earth, represented as a unit
-sphere. The rotation is specified using the roll (α, x-axis rotation), pitch (β,
-y-axis rotation), and yaw (γ, z-axis rotation), all angles are given in degrees.
-The returned point is given as a pair of lon, lat coordinates.
-"""
-function _rotate_coordinate(lon, lat, roll, pitch, yaw)
-
-    # Move everything to rads
-    α = deg2rad(roll)
-    β = deg2rad(pitch)
-    γ = deg2rad(yaw)
-
-    # Coordinates to unit sphere position
-    x = cos(deg2rad(lat)) * cos(deg2rad(lon))
-    y = cost(deg2rad(lat)) * sin(deg2rad(lon))
-    z = sin(deg2rad(lat))
-
-    # Point vector
-    p = [x, y, z]
-
-    # Rotation matrix for the roll
-    Rx = [1 0 0; 0 cos(α) -sin(α); 0 sin(α) cos(α)]
-    # Rotation matrix for the pitch
-    Ry = [cos(β) 0 sin(β); 0 1 0; -sin(β) - cos(β)]
-    # Rotation matrix for the yaw
-    Rz = [cos(γ) -sin(γ) 0; sin(γ) cos(γ) 0; 0 0 1]
-
-    # Rotation matrix (x, then y, then z)
-    R = Rz * Ry * Rx
-
-    # New point
-    X, Y, Z = R * p
-
-    # New longitude and latitude
-    projected = (atan(Y, X), asin(Z))
-    
-    # Return the new point in degrees
-    return rad2deg.(projected)
-end
+function _shift_coordinates(points::Vector{Tuple{Float64, Float64}}, angle_lon::Float64, angle_lat::Float64)
+    shifted_points = Vector{Tuple{Float64, Float64}}()
+    rotation_angle_lon = deg2rad(angle_lon)
+    rotation_angle_lat = deg2rad(angle_lat)
+
+    for (lon, lat) in points
+        # Convert longitude and latitude to radians
+        lon_rad = deg2rad(lon)
+        lat_rad = deg2rad(lat)
+
+        # Perform rotation using spherical trigonometry
+        new_lon_rad = lon_rad + rotation_angle_lon * cos(lat_rad)
+        new_lat_rad = lat_rad + rotation_angle_lat * cos(lon_rad)
+
+        # Convert back to degrees
+        new_lon = rad2deg(new_lon_rad)
+        new_lat = rad2deg(new_lat_rad)
+
+        # Ensure longitude stays within [-180, 180]
+        new_lon = mod(new_lon + 180, 360) - 180
 
-#=
-
-function generate_valid_rotation(reference, pool; roll=(-5.0, 5.0), pitch=(-5.0, 5.0), yaw=(-5.0, 5.0))
-    E, N, K = eastings(reference), northings(reference), keys(reference)
-    still_searching = true
-    Θ = (0.0, 0.0, 0.0)
-    while still_searching
-        α = rand() * (roll[2]-roll[1]) + roll[1] 
-        β = rand() * (pitch[2]-pitch[1]) + pitch[1]
-        γ = rand() * (yaw[2]-yaw[1]) + yaw[1]
-        Θ = (α, β, γ)
-        still_searching = false
-        for k in K
-            lat, lon = rotate_earth_coordinates(N[k[1]], E[k[2]], Θ...)
-            try
-                nv = pool[lon, lat]
-                if isnothing(nv)
-                    still_searching = true
-                    break
-                end
-            catch err
-                still_searching = true
-                break
-            end
-        end
+        # Ensure latitude stays within [-90, 90]
+        new_lat = max(min(new_lat, 90), -90)
+
+        # Append the rotated point to the result
+        push!(shifted_points, (new_lon, new_lat))
     end
-    return Θ
+
+    return shifted_points
 end
-=#
 
-#=
+function _rotate_coordinates(points::Vector{Tuple{Float64, Float64}}, angle::Float64)
+    # Compute the centroid of the points
+    centroid_lon = mean(map(p -> p[1], points))
+    centroid_lat = mean(map(p -> p[2], points))
 
-temp = SDMLayer(RasterData(WorldClim2, BioClim); layer=1, resolution=2.5)
-pol = getpolygon(PolygonData(NaturalEarth, Countries))["Czech Republic"]
-ref = trim(mask(temp, pol))
-lat, lon = northings(ref), eastings(ref)
+    # Convert centroid to radians
+    centroid_lon_rad = deg2rad(centroid_lon)
+    centroid_lat_rad = deg2rad(centroid_lat)
 
-# rotation, up/down, left/right
-angles = generate_valid_rotation(ref, temp)
-@info angles
+    # Convert angle to radians
+    rotation_angle = deg2rad(angle)
 
-sar = deepcopy(ref)
+    rotated_points = Vector{Tuple{Float64, Float64}}()
 
-for k in keys(ref)
-    coords = (lat[k[1]], lon[k[2]])
-    nlat, nlon = rotate_earth_coordinates(coords..., angles...)
-    sar.grid[k] = temp[nlon, nlat]
-end
+    for (lon, lat) in points
+        # Convert point to radians
+        lon_rad = deg2rad(lon)
+        lat_rad = deg2rad(lat)
+
+        # Convert to Cartesian coordinates
+        x = cos(lat_rad) * cos(lon_rad)
+        y = cos(lat_rad) * sin(lon_rad)
+        z = sin(lat_rad)
 
-replaced = quantiletransfer(sar, ref)
+        # Compute centroid in Cartesian coordinates
+        cx = cos(centroid_lat_rad) * cos(centroid_lon_rad)
+        cy = cos(centroid_lat_rad) * sin(centroid_lon_rad)
+        cz = sin(centroid_lat_rad)
 
-f = Figure()
-ax = Axis(f[1,1], aspect=DataAspect())
-ax2 = Axis(f[2,1])
-hm = heatmap!(ax, replaced, colormap=:navia)
-hist!(ax2, replaced, bins=100)
-Colorbar(f[1:2,2], hm)
-f
-=#
\ No newline at end of file
+        # Rotate around the centroid using Rodrigues' rotation formula
+        kx, ky, kz = cx, cy, cz
+        dot_product = x * kx + y * ky + z * kz
+        new_x = x * cos(rotation_angle) + (ky * z - kz * y) * sin(rotation_angle) + kx * dot_product * (1 - cos(rotation_angle))
+        new_y = y * cos(rotation_angle) + (kz * x - kx * z) * sin(rotation_angle) + ky * dot_product * (1 - cos(rotation_angle))
+        new_z = z * cos(rotation_angle) + (kx * y - ky * x) * sin(rotation_angle) + kz * dot_product * (1 - cos(rotation_angle))
+
+        # Convert back to spherical coordinates
+        new_lon_rad = atan(new_y, new_x)
+        new_lat_rad = atan(new_z, sqrt(new_x^2 + new_y^2))
+
+        # Convert back to degrees
+        new_lon = rad2deg(new_lon_rad)
+        new_lat = rad2deg(new_lat_rad)
+
+        # Ensure longitude stays within [-180, 180]
+        new_lon = mod(new_lon + 180, 360) - 180
+
+        # Ensure latitude stays within [-90, 90]
+        new_lat = max(min(new_lat, 90), -90)
+
+        # Append the rotated point to the result
+        push!(rotated_points, (new_lon, new_lat))
+    end
+
+    return rotated_points
+end

From 2385601c8a3de36ae694c4ddf3bc75426adc06c2 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Timoth=C3=A9e=20Poisot?= <timothee.poisot@umontreal.ca>
Date: Mon, 28 Jul 2025 16:10:10 -0400
Subject: [PATCH 05/19] feat(layers)?: SAR bad code to refactor

---
 SimpleSDMLayers/src/shift-and-rotate.jl | 253 +++++++++++++++++-------
 1 file changed, 184 insertions(+), 69 deletions(-)

diff --git a/SimpleSDMLayers/src/shift-and-rotate.jl b/SimpleSDMLayers/src/shift-and-rotate.jl
index c05be48ae..3c92c596a 100644
--- a/SimpleSDMLayers/src/shift-and-rotate.jl
+++ b/SimpleSDMLayers/src/shift-and-rotate.jl
@@ -1,87 +1,202 @@
-function _shift_coordinates(points::Vector{Tuple{Float64, Float64}}, angle_lon::Float64, angle_lat::Float64)
-    shifted_points = Vector{Tuple{Float64, Float64}}()
-    rotation_angle_lon = deg2rad(angle_lon)
-    rotation_angle_lat = deg2rad(angle_lat)
+#=
+using SpeciesDistributionToolkit
+const SDT = SpeciesDistributionToolkit
+using CairoMakie
+
+# Test layer
+cntrs = getpolygon(PolygonData(NaturalEarth, Countries))
+pol = cntrs["Italy"]
+P = SDMLayer(RasterData(WorldClim2, BioClim); resolution=2.5)
+L = SDMLayer(RasterData(WorldClim2, BioClim); resolution=2.5, SDT.boundingbox(pol)...)
+mask!(L, pol)
+
+# Check the data
+heatmap(L)
+
+# Step 1 - get the coordinates
+function lonlat(L::SDMLayer)
+    E, N, K = eastings(L), northings(L), keys(L)
+    xy = [(E[k[2]], N[k[1]]) for k in K]
+    return xy
+end
 
-    for (lon, lat) in points
-        # Convert longitude and latitude to radians
+# Step 2 - plot the coordinates
+scatter(lonlat(L))
+
+function _spherical_rotation(xy, degrees, axis=:z)
+    # Convert degrees to radians
+    θ = deg2rad(degrees)
+    
+    # x = cos(lat) * cos(lon)  [points toward 0°N, 0°E (Gulf of Guinea)]
+    # y = cos(lat) * sin(lon)  [points toward 0°N, 90°E (Bay of Bengal)]  
+    # z = sin(lat)             [points toward North Pole 90°N]
+    #
+    # ROTATION EFFECTS:
+    # :z rotation (around North Pole axis): 
+    #   - Keeps latitude CONSTANT, changes longitude
+    #   - Pure east-west movement along parallels
+    #   - Like spinning a globe around its axis
+    #
+    # :x rotation (around 0°N, 0°E axis):
+    #   - Changes both lat and lon
+    #   - 0°E meridian stays fixed, other meridians move
+    #   - Creates north-south arc motion for most points
+    #
+    # :y rotation (around 0°N, 90°E axis):
+    #   - Changes both lat and lon  
+    #   - 90°E meridian stays fixed, other meridians move
+    #   - Creates different north-south arc motion
+    
+    # Calculate centroid to determine rotation axis
+    n = length(xy)
+    cx = sum(p[1] for p in xy) / n  # centroid longitude
+    cy = sum(p[2] for p in xy) / n  # centroid latitude
+    
+    rotated = []
+    
+    for (lon, lat) in xy
+        # Convert spherical coordinates (lon, lat) to Cartesian (x, y, z)
+        # Assuming unit sphere (radius = 1)
         lon_rad = deg2rad(lon)
         lat_rad = deg2rad(lat)
-
-        # Perform rotation using spherical trigonometry
-        new_lon_rad = lon_rad + rotation_angle_lon * cos(lat_rad)
-        new_lat_rad = lat_rad + rotation_angle_lat * cos(lon_rad)
-
+        
+        x = cos(lat_rad) * cos(lon_rad)
+        y = cos(lat_rad) * sin(lon_rad)
+        z = sin(lat_rad)
+        
+        # Apply 3D rotation matrix based on specified axis
+        if axis == :z  # Rotation around z-axis (through poles)
+            x_rot = x * cos(θ) - y * sin(θ)
+            y_rot = x * sin(θ) + y * cos(θ)
+            z_rot = z
+        elseif axis == :y  # Rotation around y-axis
+            x_rot = x * cos(θ) + z * sin(θ)
+            y_rot = y
+            z_rot = -x * sin(θ) + z * cos(θ)
+        elseif axis == :x  # Rotation around x-axis
+            x_rot = x
+            y_rot = y * cos(θ) - z * sin(θ)
+            z_rot = y * sin(θ) + z * cos(θ)
+        else
+            error("Axis must be :x, :y, or :z")
+        end
+        
+        # Convert back to spherical coordinates
+        new_lat = asin(z_rot)
+        new_lon = atan(y_rot, x_rot)
+        
         # Convert back to degrees
-        new_lon = rad2deg(new_lon_rad)
-        new_lat = rad2deg(new_lat_rad)
-
-        # Ensure longitude stays within [-180, 180]
-        new_lon = mod(new_lon + 180, 360) - 180
-
-        # Ensure latitude stays within [-90, 90]
-        new_lat = max(min(new_lat, 90), -90)
-
-        # Append the rotated point to the result
-        push!(shifted_points, (new_lon, new_lat))
+        new_lon_deg = rad2deg(new_lon)
+        new_lat_deg = rad2deg(new_lat)
+        
+        push!(rotated, (new_lon_deg, new_lat_deg))
     end
-
-    return shifted_points
+    
+    return rotated
 end
 
-function _rotate_coordinates(points::Vector{Tuple{Float64, Float64}}, angle::Float64)
-    # Compute the centroid of the points
-    centroid_lon = mean(map(p -> p[1], points))
-    centroid_lat = mean(map(p -> p[2], points))
-
-    # Convert centroid to radians
-    centroid_lon_rad = deg2rad(centroid_lon)
-    centroid_lat_rad = deg2rad(centroid_lat)
-
-    # Convert angle to radians
-    rotation_angle = deg2rad(angle)
-
-    rotated_points = Vector{Tuple{Float64, Float64}}()
-
-    for (lon, lat) in points
-        # Convert point to radians
+# Function to rotate around an arbitrary axis passing through the centroid
+function _centroid_rotation(xy, degrees)
+    # Calculate centroid
+    n = length(xy)
+    cx = sum(p[1] for p in xy) / n  # centroid longitude
+    cy = sum(p[2] for p in xy) / n  # centroid latitude
+    
+    # Convert centroid to Cartesian coordinates for rotation axis
+    cx_rad = deg2rad(cx)
+    cy_rad = deg2rad(cy)
+    
+    # The rotation axis is the vector from Earth's center through the centroid
+    axis_x = cos(cy_rad) * cos(cx_rad)
+    axis_y = cos(cy_rad) * sin(cx_rad)
+    axis_z = sin(cy_rad)
+    
+    θ = deg2rad(degrees)
+    cos_θ = cos(θ)
+    sin_θ = sin(θ)
+    
+    rotated = []
+    
+    for (lon, lat) in xy
+        # Convert to Cartesian
         lon_rad = deg2rad(lon)
         lat_rad = deg2rad(lat)
-
-        # Convert to Cartesian coordinates
+        
         x = cos(lat_rad) * cos(lon_rad)
         y = cos(lat_rad) * sin(lon_rad)
         z = sin(lat_rad)
-
-        # Compute centroid in Cartesian coordinates
-        cx = cos(centroid_lat_rad) * cos(centroid_lon_rad)
-        cy = cos(centroid_lat_rad) * sin(centroid_lon_rad)
-        cz = sin(centroid_lat_rad)
-
-        # Rotate around the centroid using Rodrigues' rotation formula
-        kx, ky, kz = cx, cy, cz
-        dot_product = x * kx + y * ky + z * kz
-        new_x = x * cos(rotation_angle) + (ky * z - kz * y) * sin(rotation_angle) + kx * dot_product * (1 - cos(rotation_angle))
-        new_y = y * cos(rotation_angle) + (kz * x - kx * z) * sin(rotation_angle) + ky * dot_product * (1 - cos(rotation_angle))
-        new_z = z * cos(rotation_angle) + (kx * y - ky * x) * sin(rotation_angle) + kz * dot_product * (1 - cos(rotation_angle))
-
+        
+        # Rodrigues' rotation formula for rotation around arbitrary axis
+        # R = I + sin(θ)[K] + (1-cos(θ))[K]²
+        # where [K] is the skew-symmetric matrix of the rotation axis
+        
+        dot_product = axis_x * x + axis_y * y + axis_z * z
+        
+        x_rot = x * cos_θ + (axis_y * z - axis_z * y) * sin_θ + axis_x * dot_product * (1 - cos_θ)
+        y_rot = y * cos_θ + (axis_z * x - axis_x * z) * sin_θ + axis_y * dot_product * (1 - cos_θ)
+        z_rot = z * cos_θ + (axis_x * y - axis_y * x) * sin_θ + axis_z * dot_product * (1 - cos_θ)
+        
         # Convert back to spherical coordinates
-        new_lon_rad = atan(new_y, new_x)
-        new_lat_rad = atan(new_z, sqrt(new_x^2 + new_y^2))
-
+        new_lat = asin(clamp(z_rot, -1, 1))  # Clamp to handle numerical errors
+        new_lon = atan(y_rot, x_rot)
+        
         # Convert back to degrees
-        new_lon = rad2deg(new_lon_rad)
-        new_lat = rad2deg(new_lat_rad)
-
-        # Ensure longitude stays within [-180, 180]
-        new_lon = mod(new_lon + 180, 360) - 180
+        new_lon_deg = rad2deg(new_lon)
+        new_lat_deg = rad2deg(new_lat)
+        
+        push!(rotated, (new_lon_deg, new_lat_deg))
+    end
+    
+    return rotated
+end
 
-        # Ensure latitude stays within [-90, 90]
-        new_lat = max(min(new_lat, 90), -90)
+# Function for complex chained rotation attempting longitude preservation
+function _longitude_preservation(xy, latitude_shift_degrees)
+    # APPROACH 4: Multi-step rotation chain
+    # This attempts to use rotations to achieve latitude-only change
+    
+    # Calculate centroid for reference
+    n = length(xy)
+    cx = sum(p[1] for p in xy) / n
+    cy = sum(p[2] for p in xy) / n
+    
+    # Strategy: Rotate to align with a preferred meridian, shift, then rotate back
+    align_angle = -cx  # Angle to rotate to 0° longitude
+    
+    # Step 1: Rotate to align centroid with 0° meridian (Z-axis rotation)
+    aligned = rotate_sphere(xy, align_angle, :z)
+    
+    # Step 2: Apply X-axis rotation (which preserves 0° meridian)
+    shifted = rotate_sphere(aligned, latitude_shift_degrees, :y)
+    
+    # Step 3: Rotate back to original longitude position
+    return rotate_sphere(shifted, -align_angle, :z)
+end
 
-        # Append the rotated point to the result
-        push!(rotated_points, (new_lon, new_lat))
+shiftlatitude(angle) = (xy) -> _spherical_rotation(xy, angle, :z)
+shiftlongitude(angle) = (xy) -> _longitude_preservation(xy, angle)
+centroidrotation(angle) = (xy) -> _centroid_rotation(xy, angle)
+shiftandrotate(lon, lat, cntr) = (xy) -> xy |> shiftlatitude(lat) |> shiftlongitude(lon) |> centroidrotation(cntr)
+
+rangelat = LinRange(-20.0, 20.0, 40)
+rangelon = LinRange(-20.0, 20.0, 40)
+rangerot = LinRange(-20.0, 20.0, 200)
+
+r = (rand(rangelat), rand(rangelon), rand(rangerot))
+searching = true
+while searching
+    r = (rand(rangelat), rand(rangelon), rand(rangerot))
+    trf = shiftandrotate(r...)
+    u = [P[c...] for c in trf(lonlat(L))]
+    if !any(isnothing, u)
+        searching = false
     end
-
-    return rotated_points
 end
+
+@info r
+trf = shiftandrotate(r...)
+
+Y = deepcopy(L)
+SimpleSDMLayers.burnin!(Y, [P[c...] for c in trf(lonlat(L))])
+heatmap(Y)
+=#
\ No newline at end of file

From f24d210d299490e55411f4bedfd81bfe6e57d376 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Timoth=C3=A9e=20Poisot?= <timothee.poisot@umontreal.ca>
Date: Mon, 28 Jul 2025 18:32:21 -0400
Subject: [PATCH 06/19] feat(layers)?: SAR

---
 SimpleSDMLayers/src/shift-and-rotate.jl | 198 +++++++++++++-----------
 1 file changed, 109 insertions(+), 89 deletions(-)

diff --git a/SimpleSDMLayers/src/shift-and-rotate.jl b/SimpleSDMLayers/src/shift-and-rotate.jl
index 3c92c596a..9563a93d2 100644
--- a/SimpleSDMLayers/src/shift-and-rotate.jl
+++ b/SimpleSDMLayers/src/shift-and-rotate.jl
@@ -5,7 +5,7 @@ using CairoMakie
 
 # Test layer
 cntrs = getpolygon(PolygonData(NaturalEarth, Countries))
-pol = cntrs["Italy"]
+pol = cntrs["Kenya"]
 P = SDMLayer(RasterData(WorldClim2, BioClim); resolution=2.5)
 L = SDMLayer(RasterData(WorldClim2, BioClim); resolution=2.5, SDT.boundingbox(pol)...)
 mask!(L, pol)
@@ -23,164 +23,172 @@ end
 # Step 2 - plot the coordinates
 scatter(lonlat(L))
 
-function _spherical_rotation(xy, degrees, axis=:z)
+"""
+    _rotation_z(x, y, z, θ)
+
+Rotation around the polar axis -- latitude is constant but longitude changes
+"""
+function _rotation_z(x, y, z, θ)
+    x_rot = x * cos(θ) - y * sin(θ)
+    y_rot = x * sin(θ) + y * cos(θ)
+    z_rot = z
+    return (x_rot, y_rot, z_rot)
+end
+
+"""
+    _rotation_y(x, y, z, θ)
+
+Rotation along the axis ending at the Bay of Bengal
+"""
+function _rotation_y(x, y, z, θ)
+    x_rot = x * cos(θ) + z * sin(θ)
+    y_rot = y
+    z_rot = -x * sin(θ) + z * cos(θ)
+    return (x_rot, y_rot, z_rot)
+end
+
+"""
+    _rotation_x(x, y, z, θ)
+
+Rotation along the axis ending at the Bay of Guinea
+"""
+function _rotation_x(x, y, z, θ)
+    x_rot = x
+    y_rot = y * cos(θ) - z * sin(θ)
+    z_rot = y * sin(θ) + z * cos(θ)
+    return (x_rot, y_rot, z_rot)
+end
+
+"""
+    _spherical_rotation(xy, degrees, axis=1)
+
+Rotation of the Earth alongside three axes, by a given angle in degrees.
+
+Axis 1 - x rotation - rotates around the axis ending in the Bay of Guinea
+
+Axis 2 - y rotation - rotates around the axis ending in the Gulf of Bengal
+
+Axis 3 - z rotation - rotates around the axis ending at the North Pole
+"""
+function _spherical_rotation(xy, degrees, axis=1)
     # Convert degrees to radians
     θ = deg2rad(degrees)
-    
-    # x = cos(lat) * cos(lon)  [points toward 0°N, 0°E (Gulf of Guinea)]
-    # y = cos(lat) * sin(lon)  [points toward 0°N, 90°E (Bay of Bengal)]  
-    # z = sin(lat)             [points toward North Pole 90°N]
-    #
-    # ROTATION EFFECTS:
-    # :z rotation (around North Pole axis): 
-    #   - Keeps latitude CONSTANT, changes longitude
-    #   - Pure east-west movement along parallels
-    #   - Like spinning a globe around its axis
-    #
-    # :x rotation (around 0°N, 0°E axis):
-    #   - Changes both lat and lon
-    #   - 0°E meridian stays fixed, other meridians move
-    #   - Creates north-south arc motion for most points
-    #
-    # :y rotation (around 0°N, 90°E axis):
-    #   - Changes both lat and lon  
-    #   - 90°E meridian stays fixed, other meridians move
-    #   - Creates different north-south arc motion
-    
+
     # Calculate centroid to determine rotation axis
     n = length(xy)
     cx = sum(p[1] for p in xy) / n  # centroid longitude
     cy = sum(p[2] for p in xy) / n  # centroid latitude
-    
+
+    # Rotation function
+    rf = [_rotation_x, _rotation_y, _rotation_z][axis]
+
     rotated = []
-    
+
     for (lon, lat) in xy
         # Convert spherical coordinates (lon, lat) to Cartesian (x, y, z)
         # Assuming unit sphere (radius = 1)
         lon_rad = deg2rad(lon)
         lat_rad = deg2rad(lat)
-        
+
         x = cos(lat_rad) * cos(lon_rad)
         y = cos(lat_rad) * sin(lon_rad)
         z = sin(lat_rad)
-        
+
         # Apply 3D rotation matrix based on specified axis
-        if axis == :z  # Rotation around z-axis (through poles)
-            x_rot = x * cos(θ) - y * sin(θ)
-            y_rot = x * sin(θ) + y * cos(θ)
-            z_rot = z
-        elseif axis == :y  # Rotation around y-axis
-            x_rot = x * cos(θ) + z * sin(θ)
-            y_rot = y
-            z_rot = -x * sin(θ) + z * cos(θ)
-        elseif axis == :x  # Rotation around x-axis
-            x_rot = x
-            y_rot = y * cos(θ) - z * sin(θ)
-            z_rot = y * sin(θ) + z * cos(θ)
-        else
-            error("Axis must be :x, :y, or :z")
-        end
-        
+        x_rot, y_rot, z_rot = rf(x, y, z, θ)
+
         # Convert back to spherical coordinates
         new_lat = asin(z_rot)
         new_lon = atan(y_rot, x_rot)
-        
+
         # Convert back to degrees
         new_lon_deg = rad2deg(new_lon)
         new_lat_deg = rad2deg(new_lat)
-        
+
         push!(rotated, (new_lon_deg, new_lat_deg))
     end
-    
+
     return rotated
 end
 
-# Function to rotate around an arbitrary axis passing through the centroid
+"""
+    _centroid_rotation(xy, degrees)
+
+Rotates a group of points around their centroid by an angle given in degrees
+"""
 function _centroid_rotation(xy, degrees)
     # Calculate centroid
     n = length(xy)
     cx = sum(p[1] for p in xy) / n  # centroid longitude
     cy = sum(p[2] for p in xy) / n  # centroid latitude
-    
+
     # Convert centroid to Cartesian coordinates for rotation axis
     cx_rad = deg2rad(cx)
     cy_rad = deg2rad(cy)
-    
+
     # The rotation axis is the vector from Earth's center through the centroid
     axis_x = cos(cy_rad) * cos(cx_rad)
     axis_y = cos(cy_rad) * sin(cx_rad)
     axis_z = sin(cy_rad)
-    
+
     θ = deg2rad(degrees)
     cos_θ = cos(θ)
     sin_θ = sin(θ)
-    
+
     rotated = []
-    
+
     for (lon, lat) in xy
         # Convert to Cartesian
         lon_rad = deg2rad(lon)
         lat_rad = deg2rad(lat)
-        
+
         x = cos(lat_rad) * cos(lon_rad)
         y = cos(lat_rad) * sin(lon_rad)
         z = sin(lat_rad)
-        
-        # Rodrigues' rotation formula for rotation around arbitrary axis
-        # R = I + sin(θ)[K] + (1-cos(θ))[K]²
-        # where [K] is the skew-symmetric matrix of the rotation axis
-        
+
+        # Rodrigues' rotation formula
         dot_product = axis_x * x + axis_y * y + axis_z * z
-        
+
         x_rot = x * cos_θ + (axis_y * z - axis_z * y) * sin_θ + axis_x * dot_product * (1 - cos_θ)
         y_rot = y * cos_θ + (axis_z * x - axis_x * z) * sin_θ + axis_y * dot_product * (1 - cos_θ)
         z_rot = z * cos_θ + (axis_x * y - axis_y * x) * sin_θ + axis_z * dot_product * (1 - cos_θ)
-        
-        # Convert back to spherical coordinates
+
+        # Back to spherical coordinates
         new_lat = asin(clamp(z_rot, -1, 1))  # Clamp to handle numerical errors
         new_lon = atan(y_rot, x_rot)
-        
-        # Convert back to degrees
+
+        # Back to degrees
         new_lon_deg = rad2deg(new_lon)
         new_lat_deg = rad2deg(new_lat)
-        
+
         push!(rotated, (new_lon_deg, new_lat_deg))
     end
-    
+
     return rotated
 end
 
-# Function for complex chained rotation attempting longitude preservation
-function _longitude_preservation(xy, latitude_shift_degrees)
-    # APPROACH 4: Multi-step rotation chain
-    # This attempts to use rotations to achieve latitude-only change
-    
+function _rotate_latitudes(xy, degrees)
+
     # Calculate centroid for reference
     n = length(xy)
     cx = sum(p[1] for p in xy) / n
     cy = sum(p[2] for p in xy) / n
-    
-    # Strategy: Rotate to align with a preferred meridian, shift, then rotate back
+
     align_angle = -cx  # Angle to rotate to 0° longitude
-    
-    # Step 1: Rotate to align centroid with 0° meridian (Z-axis rotation)
-    aligned = rotate_sphere(xy, align_angle, :z)
-    
-    # Step 2: Apply X-axis rotation (which preserves 0° meridian)
-    shifted = rotate_sphere(aligned, latitude_shift_degrees, :y)
-    
-    # Step 3: Rotate back to original longitude position
-    return rotate_sphere(shifted, -align_angle, :z)
+
+    aligned = _spherical_rotation(xy, align_angle, 3)
+    shifted = _spherical_rotation(aligned, degrees, 2)
+    return _spherical_rotation(shifted, -align_angle, 3)
 end
 
-shiftlatitude(angle) = (xy) -> _spherical_rotation(xy, angle, :z)
-shiftlongitude(angle) = (xy) -> _longitude_preservation(xy, angle)
-centroidrotation(angle) = (xy) -> _centroid_rotation(xy, angle)
-shiftandrotate(lon, lat, cntr) = (xy) -> xy |> shiftlatitude(lat) |> shiftlongitude(lon) |> centroidrotation(cntr)
+shiftlongitudes(angle) = (xy) -> _spherical_rotation(xy, angle, 3)
+shiftlatitudes(angle) = (xy) -> _rotate_latitudes(xy, -angle)
+localrotation(angle) = (xy) -> _centroid_rotation(xy, angle)
+shiftandrotate(longitude, latitude, rotation) = (xy) -> xy |> shiftlatitudes(latitude) |> shiftlongitudes(longitude) |> localrotation(rotation)
 
-rangelat = LinRange(-20.0, 20.0, 40)
-rangelon = LinRange(-20.0, 20.0, 40)
-rangerot = LinRange(-20.0, 20.0, 200)
+rangelat = LinRange(-10.0, 10.0, 40)
+rangelon = LinRange(-10.0, 10.0, 40)
+rangerot = LinRange(-45.0, 45.0, 200)
 
 r = (rand(rangelat), rand(rangelon), rand(rangerot))
 searching = true
@@ -196,7 +204,19 @@ end
 @info r
 trf = shiftandrotate(r...)
 
+f = Figure()
+ax = Axis(f[1,1]; aspect=DataAspect())
+
+heatmap!(ax, P)
+scatter!(ax, lonlat(L))
+scatter!(ax, trf(lonlat(L)))
+xlims!(ax, extrema(first.(vcat(lonlat(L), trf(lonlat(L))))) .+ (-1, 1))
+ylims!(ax, extrema(last.(vcat(lonlat(L), trf(lonlat(L))))) .+ (-1, 1))
+
 Y = deepcopy(L)
 SimpleSDMLayers.burnin!(Y, [P[c...] for c in trf(lonlat(L))])
-heatmap(Y)
+ax2 = Axis(f[1,2]; aspect=DataAspect())
+heatmap!(ax2, Y)
+
+current_figure()
 =#
\ No newline at end of file

From 45eaca20a86b4f7ee69e3505d3205695b55f0718 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Timoth=C3=A9e=20Poisot?= <timothee.poisot@umontreal.ca>
Date: Mon, 28 Jul 2025 21:48:31 -0400
Subject: [PATCH 07/19] feat(layers): shift and rotate

---
 SimpleSDMLayers/src/SimpleSDMLayers.jl        |  7 +-
 .../rotations.jl}                             | 66 +-----------------
 .../src/shift-and-rotate/shift-and-rotate.jl  | 69 +++++++++++++++++++
 3 files changed, 76 insertions(+), 66 deletions(-)
 rename SimpleSDMLayers/src/{shift-and-rotate.jl => shift-and-rotate/rotations.jl} (70%)
 create mode 100644 SimpleSDMLayers/src/shift-and-rotate/shift-and-rotate.jl

diff --git a/SimpleSDMLayers/src/SimpleSDMLayers.jl b/SimpleSDMLayers/src/SimpleSDMLayers.jl
index 1bc8ef6ef..2cfa52aa1 100644
--- a/SimpleSDMLayers/src/SimpleSDMLayers.jl
+++ b/SimpleSDMLayers/src/SimpleSDMLayers.jl
@@ -93,6 +93,11 @@ include("quantiletransfer.jl")
 export quantiletransfer, quantiletransfer!
 
 # Shift and rotate
-include("shift-and-rotate.jl")
+include("shift-and-rotate/rotations.jl")
+include("shift-and-rotate/shift-and-rotate.jl")
+export roll, pitch, yaw
+export shiftlatitudes, shiftlongitudes, localrotation
+export shiftandrotate
+export find_rotations
 
 end # module
diff --git a/SimpleSDMLayers/src/shift-and-rotate.jl b/SimpleSDMLayers/src/shift-and-rotate/rotations.jl
similarity index 70%
rename from SimpleSDMLayers/src/shift-and-rotate.jl
rename to SimpleSDMLayers/src/shift-and-rotate/rotations.jl
index 9563a93d2..71cb5d0c5 100644
--- a/SimpleSDMLayers/src/shift-and-rotate.jl
+++ b/SimpleSDMLayers/src/shift-and-rotate/rotations.jl
@@ -1,27 +1,3 @@
-#=
-using SpeciesDistributionToolkit
-const SDT = SpeciesDistributionToolkit
-using CairoMakie
-
-# Test layer
-cntrs = getpolygon(PolygonData(NaturalEarth, Countries))
-pol = cntrs["Kenya"]
-P = SDMLayer(RasterData(WorldClim2, BioClim); resolution=2.5)
-L = SDMLayer(RasterData(WorldClim2, BioClim); resolution=2.5, SDT.boundingbox(pol)...)
-mask!(L, pol)
-
-# Check the data
-heatmap(L)
-
-# Step 1 - get the coordinates
-function lonlat(L::SDMLayer)
-    E, N, K = eastings(L), northings(L), keys(L)
-    xy = [(E[k[2]], N[k[1]]) for k in K]
-    return xy
-end
-
-# Step 2 - plot the coordinates
-scatter(lonlat(L))
 
 """
     _rotation_z(x, y, z, θ)
@@ -179,44 +155,4 @@ function _rotate_latitudes(xy, degrees)
     aligned = _spherical_rotation(xy, align_angle, 3)
     shifted = _spherical_rotation(aligned, degrees, 2)
     return _spherical_rotation(shifted, -align_angle, 3)
-end
-
-shiftlongitudes(angle) = (xy) -> _spherical_rotation(xy, angle, 3)
-shiftlatitudes(angle) = (xy) -> _rotate_latitudes(xy, -angle)
-localrotation(angle) = (xy) -> _centroid_rotation(xy, angle)
-shiftandrotate(longitude, latitude, rotation) = (xy) -> xy |> shiftlatitudes(latitude) |> shiftlongitudes(longitude) |> localrotation(rotation)
-
-rangelat = LinRange(-10.0, 10.0, 40)
-rangelon = LinRange(-10.0, 10.0, 40)
-rangerot = LinRange(-45.0, 45.0, 200)
-
-r = (rand(rangelat), rand(rangelon), rand(rangerot))
-searching = true
-while searching
-    r = (rand(rangelat), rand(rangelon), rand(rangerot))
-    trf = shiftandrotate(r...)
-    u = [P[c...] for c in trf(lonlat(L))]
-    if !any(isnothing, u)
-        searching = false
-    end
-end
-
-@info r
-trf = shiftandrotate(r...)
-
-f = Figure()
-ax = Axis(f[1,1]; aspect=DataAspect())
-
-heatmap!(ax, P)
-scatter!(ax, lonlat(L))
-scatter!(ax, trf(lonlat(L)))
-xlims!(ax, extrema(first.(vcat(lonlat(L), trf(lonlat(L))))) .+ (-1, 1))
-ylims!(ax, extrema(last.(vcat(lonlat(L), trf(lonlat(L))))) .+ (-1, 1))
-
-Y = deepcopy(L)
-SimpleSDMLayers.burnin!(Y, [P[c...] for c in trf(lonlat(L))])
-ax2 = Axis(f[1,2]; aspect=DataAspect())
-heatmap!(ax2, Y)
-
-current_figure()
-=#
\ No newline at end of file
+end
\ No newline at end of file
diff --git a/SimpleSDMLayers/src/shift-and-rotate/shift-and-rotate.jl b/SimpleSDMLayers/src/shift-and-rotate/shift-and-rotate.jl
new file mode 100644
index 000000000..2933d6b26
--- /dev/null
+++ b/SimpleSDMLayers/src/shift-and-rotate/shift-and-rotate.jl
@@ -0,0 +1,69 @@
+function lonlat(L::SDMLayer)
+    E, N, K = eastings(L), northings(L), keys(L)
+    xy = [(E[k[2]], N[k[1]]) for k in K]
+    return xy
+end
+
+
+roll(angle) = (xy) -> _spherical_rotation(xy, angle, 1)
+pitch(angle) = (xy) -> _spherical_rotation(xy, angle, 2)
+yaw(angle) = (xy) -> _spherical_rotation(xy, angle, 3)
+
+shiftlongitudes(angle) = (xy) -> _spherical_rotation(xy, angle, 3)
+shiftlatitudes(angle) = (xy) -> _rotate_latitudes(xy, -angle)
+
+localrotation(angle) = (xy) -> _centroid_rotation(xy, angle)
+
+shiftandrotate(longitude, latitude, rotation) = (xy) -> xy |> shiftlongitudes(longitude) |> shiftlatitudes(latitude) |> localrotation(rotation)
+
+
+function find_rotations(L, P; longitudes=(-10., 10.), latitudes=(-10., 10.), rotation=(-10., 10.), maxiter=10_000, steps=150)
+    iter = 1
+    rangelon = LinRange(extrema(longitudes)..., steps)
+    rangelat = LinRange(extrema(latitudes)..., steps)
+    rangerot = LinRange(extrema(rotation)..., steps)
+    r = (rand(rangelon), rand(rangelat), rand(rangerot))
+    while iter < maxiter
+        iter += 1
+        r = (rand(rangelon), rand(rangelat), rand(rangerot))
+        trf = shiftandrotate(r...)
+        u = [P[c...] for c in trf(lonlat(L))]
+        if !any(isnothing, u)
+            return r
+        end
+    end
+    return nothing
+end
+
+#=
+cntrs = getpolygon(PolygonData(NaturalEarth, Countries))
+pol = cntrs["Czech Republic"]
+P = SDMLayer(RasterData(WorldClim2, BioClim); resolution=2.5)
+L = SDMLayer(RasterData(WorldClim2, BioClim); resolution=2.5, SDT.boundingbox(pol)...)
+mask!(L, pol)
+
+# Check the data
+heatmap(L)
+
+
+r = find_rotations(L, P; latitudes=(-10., 10.), longitudes=(-20., 20.), rotation=(-1., 1.))
+@info r
+
+trf = shiftandrotate(r...)
+
+f = Figure(; size=(1000, 1000))
+ax = Axis(f[1, 1]; aspect=DataAspect())
+
+heatmap!(ax, P, colormap=:greys)
+scatter!(ax, lonlat(L), markersize=1, color=:black)
+scatter!(ax, trf(lonlat(L)), markersize=1, color=:red)
+xlims!(ax, extrema(first.(vcat(lonlat(L), trf(lonlat(L))))) .+ (-1, 1))
+ylims!(ax, extrema(last.(vcat(lonlat(L), trf(lonlat(L))))) .+ (-1, 1))
+
+Y = deepcopy(L)
+SimpleSDMLayers.burnin!(Y, [P[c...] for c in trf(lonlat(L))])
+ax2 = Axis(f[1, 2]; aspect=DataAspect())
+heatmap!(ax2, Y)
+
+current_figure()
+=#
\ No newline at end of file

From f701a7a04d5f41d1cb3521a916324d807c933e99 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Timoth=C3=A9e=20Poisot?= <timothee.poisot@umontreal.ca>
Date: Tue, 29 Jul 2025 22:16:06 -0400
Subject: [PATCH 08/19] compat(layers): LinearAlgebra

---
 SimpleSDMLayers/Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/SimpleSDMLayers/Project.toml b/SimpleSDMLayers/Project.toml
index a1133efba..97c38c1c8 100644
--- a/SimpleSDMLayers/Project.toml
+++ b/SimpleSDMLayers/Project.toml
@@ -29,7 +29,7 @@ NeutralLandscapesExtension = "NeutralLandscapes"
 ArchGDAL = "0.10"
 Clustering = "0.15"
 Distances = "0.10"
-LinearAlgebra = "1.11.0"
+LinearAlgebra = "1"
 Makie = "0.21 - 0.24"
 MultivariateStats = "0.10"
 NeutralLandscapes = "0.1"

From f308b6b3f8ff8e470e64491830f5dbc37d1dbbe6 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Timoth=C3=A9e=20Poisot?= <timothee.poisot@umontreal.ca>
Date: Tue, 29 Jul 2025 22:16:58 -0400
Subject: [PATCH 09/19] chore(layers): format

---
 SimpleSDMLayers/src/shift-and-rotate/shift-and-rotate.jl | 2 --
 1 file changed, 2 deletions(-)

diff --git a/SimpleSDMLayers/src/shift-and-rotate/shift-and-rotate.jl b/SimpleSDMLayers/src/shift-and-rotate/shift-and-rotate.jl
index 2933d6b26..47ec900cf 100644
--- a/SimpleSDMLayers/src/shift-and-rotate/shift-and-rotate.jl
+++ b/SimpleSDMLayers/src/shift-and-rotate/shift-and-rotate.jl
@@ -4,7 +4,6 @@ function lonlat(L::SDMLayer)
     return xy
 end
 
-
 roll(angle) = (xy) -> _spherical_rotation(xy, angle, 1)
 pitch(angle) = (xy) -> _spherical_rotation(xy, angle, 2)
 yaw(angle) = (xy) -> _spherical_rotation(xy, angle, 3)
@@ -16,7 +15,6 @@ localrotation(angle) = (xy) -> _centroid_rotation(xy, angle)
 
 shiftandrotate(longitude, latitude, rotation) = (xy) -> xy |> shiftlongitudes(longitude) |> shiftlatitudes(latitude) |> localrotation(rotation)
 
-
 function find_rotations(L, P; longitudes=(-10., 10.), latitudes=(-10., 10.), rotation=(-10., 10.), maxiter=10_000, steps=150)
     iter = 1
     rangelon = LinRange(extrema(longitudes)..., steps)

From e71f0ba122f7c36b16d2d2f2ff3c477e5ba3a6bd Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Timoth=C3=A9e=20Poisot?= <timothee.poisot@umontreal.ca>
Date: Tue, 29 Jul 2025 22:33:26 -0400
Subject: [PATCH 10/19] doc(layer): docstrings

---
 .../src/shift-and-rotate/shift-and-rotate.jl  | 71 +++++++++++++++++--
 1 file changed, 64 insertions(+), 7 deletions(-)

diff --git a/SimpleSDMLayers/src/shift-and-rotate/shift-and-rotate.jl b/SimpleSDMLayers/src/shift-and-rotate/shift-and-rotate.jl
index 47ec900cf..5a57f9927 100644
--- a/SimpleSDMLayers/src/shift-and-rotate/shift-and-rotate.jl
+++ b/SimpleSDMLayers/src/shift-and-rotate/shift-and-rotate.jl
@@ -1,3 +1,8 @@
+"""
+    lonlat(L::SDMLayer)
+
+Returns a vector of longitudes and latitudes for a layer.
+"""
 function lonlat(L::SDMLayer)
     E, N, K = eastings(L), northings(L), keys(L)
     xy = [(E[k[2]], N[k[1]]) for k in K]
@@ -8,22 +13,74 @@ roll(angle) = (xy) -> _spherical_rotation(xy, angle, 1)
 pitch(angle) = (xy) -> _spherical_rotation(xy, angle, 2)
 yaw(angle) = (xy) -> _spherical_rotation(xy, angle, 3)
 
+"""
+    shiftlongitudes(angle)
+
+Shifts the longitudes of a group of points by performing a rotation alonside the
+yaw axis (going through the Nort pole). This preserves the latitudes of all
+points exactly and does not deform the shape.
+"""
 shiftlongitudes(angle) = (xy) -> _spherical_rotation(xy, angle, 3)
+
+"""
+    shiftlatitudes(angle)
+    
+Returns a function to move coordinates up or down in latitudes by a given angle
+in degree. Note that this accounts for the curvature of the Earth, and therefore
+requires three distinct operations. First, the points are moved so that their
+centroid has a longitude of 0; then, the points are shifted in latitute by
+performing a rotation along the pitch axis by the desired angle; finally, the
+centroid of the points is brought back to its original longitude. For this
+reason, and because the rotation along the pitch axis accounts for the curvature
+of the Earth, this function will deform the shape, and specifically change the
+range of longitudes covered.
+"""
 shiftlatitudes(angle) = (xy) -> _rotate_latitudes(xy, -angle)
 
+"""
+    localrotation(angle)
+    
+Returns a function to create a rotation of coordinates around their centroids.
+"""
 localrotation(angle) = (xy) -> _centroid_rotation(xy, angle)
 
-shiftandrotate(longitude, latitude, rotation) = (xy) -> xy |> shiftlongitudes(longitude) |> shiftlatitudes(latitude) |> localrotation(rotation)
+"""
+    shiftandrotate(longitude::T, latitude::T, rotation::T) where T <: Number
+
+Returns a function to transform a vector of lon,lat points according to three
+operations done in order:
+
+    - a shift of the longitudes
+    - a shift of the latitudes (the order of these two operations is actually irrelevant)
+    - a local rotation around the centroid
+
+All the angles are given in degrees. The output of this function is a function
+that takes a vector of coordinates to transform. The transformations do account
+for the curvature of the Earth. For this reason, rotations and changes in the
+latitudes will deform the points, but shift in the longitudes will not.
+"""
+function shiftandrotate(longitude::T, latitude::T, rotation::T) where T <: Number
+    return shiftlongitudes(longitude) ∘ shiftlatitudes(latitude) ∘ localrotation(rotation)
+end
+
+"""
+    findrotation(L, P; longitudes=-10:0.1:10, latitudes=-10:0.1:10, rotations=-10:0.1:10, maxiter=10_000)
+
+Find a possible rotation for the `shiftandrotate` function, by attempting to
+move a target layer `L` until all of the shifted and rotated coordinates are
+valid coordinates in the layer `P`. The range of angles to explore is given as
+keywords, and the function accepts a `maxiter` argument after which, if no
+rotation is found, it returns `nothing`.
 
-function find_rotations(L, P; longitudes=(-10., 10.), latitudes=(-10., 10.), rotation=(-10., 10.), maxiter=10_000, steps=150)
+Note that it is almost always a valid strategy to look for shifts and rotations
+on a raster at a coarser resolution.
+"""
+function findrotation(L::SDMLayer, P::SDMLayer; longitudes=-10:0.1:10, latitudes=-10:0.1:10, rotations=-10:0.1:10, maxiter=10_000)
     iter = 1
-    rangelon = LinRange(extrema(longitudes)..., steps)
-    rangelat = LinRange(extrema(latitudes)..., steps)
-    rangerot = LinRange(extrema(rotation)..., steps)
-    r = (rand(rangelon), rand(rangelat), rand(rangerot))
+    r = (rand(longitudes), rand(latitudes), rand(rotations))
     while iter < maxiter
         iter += 1
-        r = (rand(rangelon), rand(rangelat), rand(rangerot))
+        r = (rand(longitudes), rand(latitudes), rand(rotations))
         trf = shiftandrotate(r...)
         u = [P[c...] for c in trf(lonlat(L))]
         if !any(isnothing, u)

From 3a45d45a03ff04616baa4768c54b8a3098150404 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Timoth=C3=A9e=20Poisot?= <timothee.poisot@umontreal.ca>
Date: Tue, 29 Jul 2025 22:37:22 -0400
Subject: [PATCH 11/19] feat(layers): lonlat + additional docstrings

---
 SimpleSDMLayers/src/SimpleSDMLayers.jl                   | 3 ++-
 SimpleSDMLayers/src/shift-and-rotate/shift-and-rotate.jl | 8 ++++++--
 2 files changed, 8 insertions(+), 3 deletions(-)

diff --git a/SimpleSDMLayers/src/SimpleSDMLayers.jl b/SimpleSDMLayers/src/SimpleSDMLayers.jl
index 2cfa52aa1..7e400d547 100644
--- a/SimpleSDMLayers/src/SimpleSDMLayers.jl
+++ b/SimpleSDMLayers/src/SimpleSDMLayers.jl
@@ -95,9 +95,10 @@ export quantiletransfer, quantiletransfer!
 # Shift and rotate
 include("shift-and-rotate/rotations.jl")
 include("shift-and-rotate/shift-and-rotate.jl")
+export lonlat
 export roll, pitch, yaw
 export shiftlatitudes, shiftlongitudes, localrotation
 export shiftandrotate
-export find_rotations
+export findrotation
 
 end # module
diff --git a/SimpleSDMLayers/src/shift-and-rotate/shift-and-rotate.jl b/SimpleSDMLayers/src/shift-and-rotate/shift-and-rotate.jl
index 5a57f9927..64caa287b 100644
--- a/SimpleSDMLayers/src/shift-and-rotate/shift-and-rotate.jl
+++ b/SimpleSDMLayers/src/shift-and-rotate/shift-and-rotate.jl
@@ -1,11 +1,15 @@
 """
     lonlat(L::SDMLayer)
 
-Returns a vector of longitudes and latitudes for a layer.
+Returns a vector of longitudes and latitudes for a layer. This will handle the
+CRS of the layer correctly. Note that only the positions that are valued (_i.e._
+using `keys`) will be returned. The values are given in an order suitable for
+use in `burnin!`.
 """
 function lonlat(L::SDMLayer)
     E, N, K = eastings(L), northings(L), keys(L)
-    xy = [(E[k[2]], N[k[1]]) for k in K]
+    prj = SimpleSDMLayers.Proj.Transformation(L.crs, "+proj=longlat +datum=WGS84 +no_defs"; always_xy=true)
+    xy = [prj(E[k[2]], N[k[1]]) for k in K]    
     return xy
 end
 

From 38ac8fc8cbee8915602d42ba0538a3ced89e2e5c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Timoth=C3=A9e=20Poisot?= <timothee.poisot@umontreal.ca>
Date: Tue, 29 Jul 2025 22:40:38 -0400
Subject: [PATCH 12/19] feat(layers): use all values for quantile transfer

---
 SimpleSDMLayers/src/quantiletransfer.jl | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/SimpleSDMLayers/src/quantiletransfer.jl b/SimpleSDMLayers/src/quantiletransfer.jl
index 7c27777ed..5b019555a 100644
--- a/SimpleSDMLayers/src/quantiletransfer.jl
+++ b/SimpleSDMLayers/src/quantiletransfer.jl
@@ -1,13 +1,13 @@
 """
-    quantiletransfer!(target::SDMLayer{T}, reference::SDMLayer; n::Int=10) where {T <: Real}
+    quantiletransfer!(target::SDMLayer{T}, reference::SDMLayer) where {T <: Real}
 
 Replaces the values in the `target` layer so that they follow the distribution
 of the values in the `reference` layer. This works by (i) identifying the
 quantile of each cell in the target layer, then (ii) replacing this with the
 value for the same quantile in the reference layer.
 """
-function quantiletransfer!(target::SDMLayer{T}, reference::SDMLayer; n::Int=10) where {T <: Real}
-    Qt = quantize(target, round(Int, n * sqrt(count(reference))))
+function quantiletransfer!(target::SDMLayer{T}, reference::SDMLayer{T}) where {T <: Real}
+    Qt = quantize(target)
     target.grid = quantile(reference, Qt.grid)
     return target
 end
@@ -17,7 +17,7 @@ end
 
 Non-mutating version of `quantiletransfer!`
 """
-function quantiletransfer(target, reference)
+function quantiletransfer(target::SDMLayer{T}, reference::SDMLayer{T}) where {T <: Real}
     t = deepcopy(target)
     quantiletransfer!(t, reference)
     return t

From 7aa4fb01d9ca68cd34d53f29d3cd9df7ad30f499 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Timoth=C3=A9e=20Poisot?= <timothee.poisot@umontreal.ca>
Date: Wed, 30 Jul 2025 09:52:49 -0400
Subject: [PATCH 13/19] doc?: add shift and rotate vignette

---
 docs/src/.vitepress/config.mts         |  1 +
 docs/src/tutorials/shift-and-rotate.jl | 52 ++++++++++++++++++++++++++
 2 files changed, 53 insertions(+)
 create mode 100644 docs/src/tutorials/shift-and-rotate.jl

diff --git a/docs/src/.vitepress/config.mts b/docs/src/.vitepress/config.mts
index 27e2c67d9..a89f24303 100644
--- a/docs/src/.vitepress/config.mts
+++ b/docs/src/.vitepress/config.mts
@@ -179,6 +179,7 @@ export default defineConfig({
                         { text: "Building a species distribution model", link: "/tutorials/sdm-training/" },
                         { text: "Ensemble models", link: "/tutorials/sdm-ensembles/" },
                         { text: "SDM with AdaBoost", link: "/tutorials/sdm-boosting/" },
+                        { text: "Shift and rotate", link: "/tutorials/shift-and-rotate/" },
                     ]
                 }
             ],
diff --git a/docs/src/tutorials/shift-and-rotate.jl b/docs/src/tutorials/shift-and-rotate.jl
new file mode 100644
index 000000000..e9a2e217f
--- /dev/null
+++ b/docs/src/tutorials/shift-and-rotate.jl
@@ -0,0 +1,52 @@
+# # Shift and rotate
+
+# In this tutorial, we will apply the [Shift and
+# Rotate](https://onlinelibrary.wiley.com/doi/abs/10.1111/2041-210X.14443)
+# technique to generate realistic maps of predictors to provide a null sample
+# for the performance of an SDM. We will also look at the quantile mapping
+# function, which ensures that the null sample has the exact same distribution
+# of values compared to the original data.
+
+using SpeciesDistributionToolkit
+const SDT = SpeciesDistributionToolkit
+using CairoMakie
+CairoMakie.activate!(; type = "png", px_per_unit = 2) #hide
+import Random #hide
+Random.seed!(123451234123121); #hide
+
+# We will get some occurrences from the butterfly *Aglais caschmirensis* in Pakistan:
+
+pol = getpolygon(PolygonData(NaturalEarth, Countries))["Pakistan"]
+presences = GBIF.download("10.15468/dl.emv5tj");
+extent = SDT.boundingbox(pol)
+
+# And then get some bioclimatic variables:
+
+provider = RasterData(CHELSA2, BioClim)
+L = SDMLayer{Float32}[
+    SDMLayer(
+        provider;
+        layer = x,
+        extent...,
+    ) for x in 1:19
+];
+mask!(L, pol)
+
+# As for the previous tutorials, we will generate some background points:
+
+presencelayer = mask(first(L), Occurrences(mask(presences, pol)))
+background = pseudoabsencemask(DistanceToEvent, presencelayer)
+bgpoints = backgroundpoints(nodata(background, d -> d < 20), 2sum(presencelayer))
+
+# And finally train an ensemble model:
+
+sdm = Bagging(SDM(RawData, DecisionTree, L, presencelayer, bgpoints), 50)
+bagfeatures!(sdm)
+
+# Let's look at the first prediction of this model:
+
+train!(sdm)
+
+heatmap(predict(sdm, L))
+lines!(pol)
+current_figure()
\ No newline at end of file

From 2e3759f55eea22af5b744019f707982b5d3165c5 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Timoth=C3=A9e=20Poisot?= <timothee.poisot@umontreal.ca>
Date: Mon, 3 Nov 2025 09:44:43 -0500
Subject: [PATCH 14/19] doc: virtual species parameter fix

---
 docs/src/tutorials/virtual-species.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/src/tutorials/virtual-species.jl b/docs/src/tutorials/virtual-species.jl
index f1964ff80..699be6de0 100644
--- a/docs/src/tutorials/virtual-species.jl
+++ b/docs/src/tutorials/virtual-species.jl
@@ -58,7 +58,7 @@ function virtualspecies(L::Vector{<:SDMLayer}; prevalence::Float64 = 0.25)
     while invalid
         centers = [rand(values(l)) for l in L]
         shapes = randn(length(L))
-        predictors = [logistic(centers[i], shapes[i]) for i in eachindex(L)]
+        predictors = [logistic(shapes[i], centers[i]) for i in eachindex(L)]
         predictions = [predictors[i].(L[i]) for i in eachindex(L)]
         rescale!.(predictions)
         global prediction = prod(predictions)

From 700ed6e5de9087c35af420bbd46bdad775ae9eb2 Mon Sep 17 00:00:00 2001
From: tpoisot <timothee.poisot@umontreal.ca>
Date: Tue, 4 Nov 2025 13:21:00 -0500
Subject: [PATCH 15/19] doc: shift

---
 docs/src/tutorials/shift-and-rotate.jl | 68 +++++++++++++++++++++++---
 1 file changed, 60 insertions(+), 8 deletions(-)

diff --git a/docs/src/tutorials/shift-and-rotate.jl b/docs/src/tutorials/shift-and-rotate.jl
index e9a2e217f..0593d83cd 100644
--- a/docs/src/tutorials/shift-and-rotate.jl
+++ b/docs/src/tutorials/shift-and-rotate.jl
@@ -38,15 +38,67 @@ presencelayer = mask(first(L), Occurrences(mask(presences, pol)))
 background = pseudoabsencemask(DistanceToEvent, presencelayer)
 bgpoints = backgroundpoints(nodata(background, d -> d < 20), 2sum(presencelayer))
 
-# And finally train an ensemble model:
+# And finally train a model:
 
-sdm = Bagging(SDM(RawData, DecisionTree, L, presencelayer, bgpoints), 50)
-bagfeatures!(sdm)
+sdm = SDM(RawData, Logistic, L, presencelayer, bgpoints)
+variables!(sdm, ForwardSelection)
 
-# Let's look at the first prediction of this model:
+# Now, we can make the prediction based on this model
 
-train!(sdm)
+# fig-shift-heatmap
+fig = Figure(; size = (900, 700))
+panel = Axis(
+    fig[1, 1];
+    xlabel = "Easting",
+    ylabel = "Northing",
+    aspect = DataAspect(),
+)
+hidedecorations!(panel)
+hidespines!(panel)
+heatmap!(
+    panel,
+    predict(sdm, L),
+    colormap=[:white, :purple]
+)
+lines!(panel, pol, color=:black)
+current_figure() #hide
 
-heatmap(predict(sdm, L))
-lines!(pol)
-current_figure()
\ No newline at end of file
+# This represents the prediction on the original data. Now, we will identify a
+# rotation under which we can generate a different landscape. Before we do this,
+# we need a pool of layers to get the potential values from, so we will use a
+# version with a padding:
+
+poolextent = SDT.boundingbox(pol; padding=10.0)
+P = SDMLayer{Float32}[
+    SDMLayer(
+        provider;
+        layer = x,
+        poolextent...,
+    ) for x in 1:19
+];
+
+# We now find a rotation:
+
+θ = findrotation(first(L), first(P))
+
+# The rotation is a shift in latitude and longitude, and then a rotation around
+# the axis.
+
+# It can be turned into a function to generate a new layer. 
+
+sar = shiftandrotate(θ...)
+
+# We can visualize the output of this transformation:
+
+# fig-shift-transform
+f = Figure()
+ax = Axis(f[1, 1]; aspect=DataAspect())
+
+heatmap!(ax, P, colormap=:greys)
+scatter!(ax, lonlat(L), markersize=1, color=:black)
+scatter!(ax, trf(lonlat(L)), markersize=1, color=:red)
+xlims!(ax, extrema(first.(vcat(lonlat(L), trf(lonlat(L))))) .+ (-1, 1))
+ylims!(ax, extrema(last.(vcat(lonlat(L), trf(lonlat(L))))) .+ (-1, 1))
+current_figure() #hide
+
+# We can then put the 
\ No newline at end of file

From 87e08a9ff41b61830d81aee61ed55a73ae03484d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Timoth=C3=A9e=20Poisot?= <timothee.poisot@umontreal.ca>
Date: Tue, 4 Nov 2025 13:36:54 -0500
Subject: [PATCH 16/19] doc: sar

---
 docs/src/tutorials/shift-and-rotate.jl | 2 --
 1 file changed, 2 deletions(-)

diff --git a/docs/src/tutorials/shift-and-rotate.jl b/docs/src/tutorials/shift-and-rotate.jl
index 0593d83cd..4193db7be 100644
--- a/docs/src/tutorials/shift-and-rotate.jl
+++ b/docs/src/tutorials/shift-and-rotate.jl
@@ -100,5 +100,3 @@ scatter!(ax, trf(lonlat(L)), markersize=1, color=:red)
 xlims!(ax, extrema(first.(vcat(lonlat(L), trf(lonlat(L))))) .+ (-1, 1))
 ylims!(ax, extrema(last.(vcat(lonlat(L), trf(lonlat(L))))) .+ (-1, 1))
 current_figure() #hide
-
-# We can then put the 
\ No newline at end of file

From 1c1af5bfdb7cf50743a102fa82b590c9e4c3a405 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Timoth=C3=A9e=20Poisot?= <timothee.poisot@umontreal.ca>
Date: Tue, 4 Nov 2025 13:48:11 -0500
Subject: [PATCH 17/19] sar vignette

---
 SimpleSDMLayers/src/SimpleSDMLayers.jl        |  6 +-
 .../src/shift-and-rotate/shift-and-rotate.jl  | 61 +++++++++----------
 docs/src/tutorials/shift-and-rotate.jl        | 54 +++++++++++++++-
 3 files changed, 83 insertions(+), 38 deletions(-)

diff --git a/SimpleSDMLayers/src/SimpleSDMLayers.jl b/SimpleSDMLayers/src/SimpleSDMLayers.jl
index 7e400d547..2ed2292b0 100644
--- a/SimpleSDMLayers/src/SimpleSDMLayers.jl
+++ b/SimpleSDMLayers/src/SimpleSDMLayers.jl
@@ -96,9 +96,9 @@ export quantiletransfer, quantiletransfer!
 include("shift-and-rotate/rotations.jl")
 include("shift-and-rotate/shift-and-rotate.jl")
 export lonlat
-export roll, pitch, yaw
-export shiftlatitudes, shiftlongitudes, localrotation
+# export roll, pitch, yaw
+# export shiftlatitudes, shiftlongitudes, localrotation
+export findrotation, rotator
 export shiftandrotate
-export findrotation
 
 end # module
diff --git a/SimpleSDMLayers/src/shift-and-rotate/shift-and-rotate.jl b/SimpleSDMLayers/src/shift-and-rotate/shift-and-rotate.jl
index 64caa287b..1aecc4cb8 100644
--- a/SimpleSDMLayers/src/shift-and-rotate/shift-and-rotate.jl
+++ b/SimpleSDMLayers/src/shift-and-rotate/shift-and-rotate.jl
@@ -49,7 +49,7 @@ Returns a function to create a rotation of coordinates around their centroids.
 localrotation(angle) = (xy) -> _centroid_rotation(xy, angle)
 
 """
-    shiftandrotate(longitude::T, latitude::T, rotation::T) where T <: Number
+    rotator(longitude::T, latitude::T, rotation::T) where T <: Number
 
 Returns a function to transform a vector of lon,lat points according to three
 operations done in order:
@@ -63,7 +63,7 @@ that takes a vector of coordinates to transform. The transformations do account
 for the curvature of the Earth. For this reason, rotations and changes in the
 latitudes will deform the points, but shift in the longitudes will not.
 """
-function shiftandrotate(longitude::T, latitude::T, rotation::T) where T <: Number
+function rotator(longitude::T, latitude::T, rotation::T) where T <: Number
     return shiftlongitudes(longitude) ∘ shiftlatitudes(latitude) ∘ localrotation(rotation)
 end
 
@@ -76,6 +76,9 @@ valid coordinates in the layer `P`. The range of angles to explore is given as
 keywords, and the function accepts a `maxiter` argument after which, if no
 rotation is found, it returns `nothing`.
 
+The output of this method is either `nothing` (no valid rotation found), or a
+rotator function (see the documentation for `rotator`).
+
 Note that it is almost always a valid strategy to look for shifts and rotations
 on a raster at a coarser resolution.
 """
@@ -85,44 +88,36 @@ function findrotation(L::SDMLayer, P::SDMLayer; longitudes=-10:0.1:10, latitudes
     while iter < maxiter
         iter += 1
         r = (rand(longitudes), rand(latitudes), rand(rotations))
-        trf = shiftandrotate(r...)
+        trf = rotator(r...)
         u = [P[c...] for c in trf(lonlat(L))]
         if !any(isnothing, u)
-            return r
+            trf
         end
     end
     return nothing
 end
 
-#=
-cntrs = getpolygon(PolygonData(NaturalEarth, Countries))
-pol = cntrs["Czech Republic"]
-P = SDMLayer(RasterData(WorldClim2, BioClim); resolution=2.5)
-L = SDMLayer(RasterData(WorldClim2, BioClim); resolution=2.5, SDT.boundingbox(pol)...)
-mask!(L, pol)
-
-# Check the data
-heatmap(L)
-
-
-r = find_rotations(L, P; latitudes=(-10., 10.), longitudes=(-20., 20.), rotation=(-1., 1.))
-@info r
-
-trf = shiftandrotate(r...)
-
-f = Figure(; size=(1000, 1000))
-ax = Axis(f[1, 1]; aspect=DataAspect())
+"""
+    shiftandrotate!(L::SDMLayer, P::SDMLayer, rotator)
 
-heatmap!(ax, P, colormap=:greys)
-scatter!(ax, lonlat(L), markersize=1, color=:black)
-scatter!(ax, trf(lonlat(L)), markersize=1, color=:red)
-xlims!(ax, extrema(first.(vcat(lonlat(L), trf(lonlat(L))))) .+ (-1, 1))
-ylims!(ax, extrema(last.(vcat(lonlat(L), trf(lonlat(L))))) .+ (-1, 1))
+Returns the output of shifting and rotating the layer `L` on the layer `P`,
+according to a `rotator` function, obtained with for example `findrotation`.
+This function overwrites `L`.
+"""
+function shiftandrotate!(L::SDMLayer, P::SDMLayer, rotator)
+    SimpleSDMLayers.burnin!(L, [P[c...] for c in rotator(lonlat(L))])
+    return L
+end
 
-Y = deepcopy(L)
-SimpleSDMLayers.burnin!(Y, [P[c...] for c in trf(lonlat(L))])
-ax2 = Axis(f[1, 2]; aspect=DataAspect())
-heatmap!(ax2, Y)
+"""
+    shiftandrotate(L::SDMLayer, P::SDMLayer, rotator)
 
-current_figure()
-=#
\ No newline at end of file
+Returns the output of shifting and rotating the layer `L` on the layer `P`,
+according to a `rotator` function, obtained with for example `findrotation`.
+This function generates a copy of `L`.
+"""
+function shiftandrotate(L::SDMLayer, P::SDMLayer, rotator)
+    Y = similar(L)
+    shiftandrotate!(Y, P, rotator)
+    return Y
+end
\ No newline at end of file
diff --git a/docs/src/tutorials/shift-and-rotate.jl b/docs/src/tutorials/shift-and-rotate.jl
index 4193db7be..0f0208e54 100644
--- a/docs/src/tutorials/shift-and-rotate.jl
+++ b/docs/src/tutorials/shift-and-rotate.jl
@@ -84,9 +84,13 @@ P = SDMLayer{Float32}[
 # The rotation is a shift in latitude and longitude, and then a rotation around
 # the axis.
 
-# It can be turned into a function to generate a new layer. 
+# It can be turned into a function to generate the rotation:
 
-sar = shiftandrotate(θ...)
+rotator(θ...)
+
+# We can, for example, shift and rotate the first layer:
+
+Y = shiftandrotate(first(L), first(P), rotator(θ...))
 
 # We can visualize the output of this transformation:
 
@@ -100,3 +104,49 @@ scatter!(ax, trf(lonlat(L)), markersize=1, color=:red)
 xlims!(ax, extrema(first.(vcat(lonlat(L), trf(lonlat(L))))) .+ (-1, 1))
 ylims!(ax, extrema(last.(vcat(lonlat(L), trf(lonlat(L))))) .+ (-1, 1))
 current_figure() #hide
+
+# Now, also ensure that although the spatial structure of this layer has
+# changed, it should have the same distribution of values as the original layer.
+# This can be done with the quantile transfer function:
+
+quantiletransfer!(Y, L)
+
+# We can check the extrema of both layers:
+
+extrema(Y)
+
+# and the original layer
+
+extrema(L)
+
+# All things together, we can now create a series of shifted and rotated layers
+# with the original distribution of values:
+
+Z = [
+    quantiletransfer!(
+        shiftandrotate(L[i], P[i], rotator(θ...)), 
+        L[i])
+    for i in eachindex(L)
+]
+
+# We can apply the previously trained SDM on the new layers
+
+# fig-shift-heatmap2
+fig = Figure(; size = (900, 700))
+panel = Axis(
+    fig[1, 1];
+    xlabel = "Easting",
+    ylabel = "Northing",
+    aspect = DataAspect(),
+)
+hidedecorations!(panel)
+hidespines!(panel)
+heatmap!(
+    panel,
+    predict(sdm, Z),
+    colormap=[:white, :purple]
+)
+lines!(panel, pol, color=:black)
+current_figure() #hide
+
+# We can, of course, also decide to train a different model based on these data, etc.
\ No newline at end of file

From 9235b87d7e31940fcfc76f698149edf28980c2f6 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Timoth=C3=A9e=20Poisot?= <timothee.poisot@umontreal.ca>
Date: Tue, 4 Nov 2025 14:48:13 -0500
Subject: [PATCH 18/19] doc: what is wrong with gbif today?

---
 docs/src/howto/get-boundingbox.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/docs/src/howto/get-boundingbox.jl b/docs/src/howto/get-boundingbox.jl
index 45d1d8e3e..fc367ba4c 100644
--- a/docs/src/howto/get-boundingbox.jl
+++ b/docs/src/howto/get-boundingbox.jl
@@ -7,7 +7,7 @@ CairoMakie.activate!(; type = "png", px_per_unit = 2) #hide
 # The `boundingbox` method accepts most types that `SpeciesDistributionToolkit`
 # knows about, and returns a tuple:
 
-occ = occurrences("hasCoordinate" => true, "country" => "CH", "limit" => 300)
+occ = occurrences("hasCoordinate" => true, "country" => "AT", "limit" => 300)
 SpeciesDistributionToolkit.boundingbox(occ)
 
 # We can also specify a padding (in degrees) that is added to the bounding box:
@@ -44,8 +44,8 @@ current_figure() #hide
 # The same method also applies to polygons:
 
 cnt = PolygonData(OpenStreetMap, Countries)
-CHE = getpolygon(cnt, country="Switzerland")
-SpeciesDistributionToolkit.boundingbox(CHE; padding = 0.5)
+pol = getpolygon(cnt, country="Austria")
+SpeciesDistributionToolkit.boundingbox(pol; padding = 0.5)
 
 # ```@meta
 # CollapsedDocStrings = true

From 61ddd166e0a096970eb6a342a47797b86de0fd05 Mon Sep 17 00:00:00 2001
From: tpoisot <timothee.poisot@umontreal.ca>
Date: Tue, 4 Nov 2025 15:36:55 -0500
Subject: [PATCH 19/19] perf(layers): SAR pre-calc of latlon

---
 SimpleSDMLayers/src/shift-and-rotate/shift-and-rotate.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/SimpleSDMLayers/src/shift-and-rotate/shift-and-rotate.jl b/SimpleSDMLayers/src/shift-and-rotate/shift-and-rotate.jl
index 1aecc4cb8..825930e40 100644
--- a/SimpleSDMLayers/src/shift-and-rotate/shift-and-rotate.jl
+++ b/SimpleSDMLayers/src/shift-and-rotate/shift-and-rotate.jl
@@ -85,11 +85,12 @@ on a raster at a coarser resolution.
 function findrotation(L::SDMLayer, P::SDMLayer; longitudes=-10:0.1:10, latitudes=-10:0.1:10, rotations=-10:0.1:10, maxiter=10_000)
     iter = 1
     r = (rand(longitudes), rand(latitudes), rand(rotations))
+    ll = lonlat(L)
     while iter < maxiter
         iter += 1
         r = (rand(longitudes), rand(latitudes), rand(rotations))
         trf = rotator(r...)
-        u = [P[c...] for c in trf(lonlat(L))]
+        u = [P[c...] for c in trf(ll)]
         if !any(isnothing, u)
             trf
         end