diff --git a/SimpleSDMLayers/CHANGELOG.md b/SimpleSDMLayers/CHANGELOG.md
index 9730ff0f9..b8bed4535 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.2`
 
 - **added** a `Makie` compatibility entry for `0.24`
diff --git a/SimpleSDMLayers/Project.toml b/SimpleSDMLayers/Project.toml
index b86cc9f48..97c38c1c8 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.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"
diff --git a/SimpleSDMLayers/src/SimpleSDMLayers.jl b/SimpleSDMLayers/src/SimpleSDMLayers.jl
index de9647e5e..2ed2292b0 100644
--- a/SimpleSDMLayers/src/SimpleSDMLayers.jl
+++ b/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 findrotation, rotator
+export shiftandrotate
+
 end # module
diff --git a/SimpleSDMLayers/src/quantiletransfer.jl b/SimpleSDMLayers/src/quantiletransfer.jl
new file mode 100644
index 000000000..5b019555a
--- /dev/null
+++ b/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
\ No newline at end of file
diff --git a/SimpleSDMLayers/src/shift-and-rotate/rotations.jl b/SimpleSDMLayers/src/shift-and-rotate/rotations.jl
new file mode 100644
index 000000000..71cb5d0c5
--- /dev/null
+++ b/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
\ 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..825930e40
--- /dev/null
+++ b/SimpleSDMLayers/src/shift-and-rotate/shift-and-rotate.jl
@@ -0,0 +1,124 @@
+"""
+    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)
+
+"""
+    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:
+
+    - 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 rotator(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`.
+
+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.
+"""
+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(ll)]
+        if !any(isnothing, u)
+            trf
+        end
+    end
+    return nothing
+end
+
+"""
+    shiftandrotate!(L::SDMLayer, P::SDMLayer, rotator)
+
+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
+
+"""
+    shiftandrotate(L::SDMLayer, P::SDMLayer, rotator)
+
+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/.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/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
diff --git a/docs/src/tutorials/shift-and-rotate.jl b/docs/src/tutorials/shift-and-rotate.jl
new file mode 100644
index 000000000..0f0208e54
--- /dev/null
+++ b/docs/src/tutorials/shift-and-rotate.jl
@@ -0,0 +1,152 @@
+# # 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 a model:
+
+sdm = SDM(RawData, Logistic, L, presencelayer, bgpoints)
+variables!(sdm, ForwardSelection)
+
+# Now, we can make the prediction based on this model
+
+# 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
+
+# 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 the rotation:
+
+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:
+
+# 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
+
+# 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
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)