Shift and rotate

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.2`
 
 - **added** a `Makie` compatibility entry for `0.24`

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.2"
+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"
 Makie = "0.21 - 0.24"
 MultivariateStats = "0.10"
 NeutralLandscapes = "0.1"

SimpleSDMLayers/src/SimpleSDMLayers.jl

@@ -9,6 +9,8 @@ rasters.
 """
 module SimpleSDMLayers
 
+using LinearAlgebra
+using Statistics
 using TestItems
 
 import ArchGDAL
@@ -86,4 +88,17 @@ export cellarea
 include("reclassify.jl")
 export reclassify
 
+# Quantile transfer
+include("quantiletransfer.jl")
+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 findrotation
+
 end # module

SimpleSDMLayers/src/quantiletransfer.jl

@@ -0,0 +1,24 @@
+"""
+    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{T}) where {T <: Real}
+    Qt = quantize(target)
+    target.grid = quantile(reference, Qt.grid)
+    return target
+end
+
+"""
+    quantiletransfer(target, reference)
+
+Non-mutating version of `quantiletransfer!`
+"""
+function quantiletransfer(target::SDMLayer{T}, reference::SDMLayer{T}) where {T <: Real}
+    t = deepcopy(target)
+    quantiletransfer!(t, reference)
+    return t
+end

SimpleSDMLayers/src/shift-and-rotate/rotations.jl

@@ -0,0 +1,158 @@
+
+"""
+    _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)
+
+    # 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
+        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
+
+"""
+    _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
+        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_θ)
+
+        # Back to spherical coordinates
+        new_lat = asin(clamp(z_rot, -1, 1))  # Clamp to handle numerical errors
+        new_lon = atan(y_rot, x_rot)
+
+        # 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 _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
+
+    align_angle = -cx  # Angle to rotate to 0° longitude
+
+    aligned = _spherical_rotation(xy, align_angle, 3)
+    shifted = _spherical_rotation(aligned, degrees, 2)
+    return _spherical_rotation(shifted, -align_angle, 3)
+end

SimpleSDMLayers/src/shift-and-rotate/shift-and-rotate.jl

@@ -0,0 +1,128 @@
+"""
+    lonlat(L::SDMLayer)
+
+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)
+    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
+
+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::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`.
+
+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
+    r = (rand(longitudes), rand(latitudes), rand(rotations))
+    while iter < maxiter
+        iter += 1
+        r = (rand(longitudes), rand(latitudes), rand(rotations))
+        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()
+=#

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/" },
                     ]
                 }
             ],