Interpolating to a new projection

using SpeciesDistributionToolkit
using CairoMakie

The interpolate method can be used to project data into another coordinate system. For example, we can get bird richness data for metropolitan France and Corsica, in the Eckert IV projection:

spatial_extent = (; left = -4.87, right = 9.63, bottom = 41.31, top = 51.14)
dataprovider = RasterData(BiodiversityMapping, BirdRichness)
layer = SDMLayer(dataprovider; layer = "Birds", spatial_extent...)

🗺️  A 104 × 123 layer with 9579 UInt16 cells
   Spatial Reference System: +proj=eck4 +lon_0=0 +x_0=0 +y_0=0 + ... _defs

The projection of a layer can be seen with:

projection(layer)

Spatial Reference System: +proj=eck4 +lon_0=0 +x_0=0 +y_0=0 + ... _defs

Re-projecting other objects

Objects that are not layers can be projected for data visualization with reproject – there is a manual entry for this use-case.

We can check out the original data:

Code for the figure

fig, ax, hm = heatmap(
    layer;
    colormap = :navia,
    figure = (; size = (800, 400)),
    axis = (; aspect = DataAspect()),
)
Colorbar(fig[:, end + 1], hm)

And project them to the more locally appropriate EPSG:27574:

ws = interpolate(layer; dest = "EPSG:27574")

🗺️  A 104 × 123 layer with 8678 Float32 cells
   Spatial Reference System: +proj=lcc +lat_1=42.165 +lat_0=42.1 ... _defs

Note that we are specifying the projection as an EPSG code, but the package is agnostic as to the way to specify it. Using a WKT or a PROJ4 string would work just as well. This is because, internally, the projection string will be parsed and converted to an appropriate representation. This is generally the PROJ4 string for coordinates projection, and the WKT for export.

By default, this produces a layer with the same dimension as the input, and uses bilinear interpolation:

Code for the figure

fig, ax, hm = heatmap(
    ws;
    colormap = :navia,
    figure = (; size = (800, 400)),
    axis = (; aspect = DataAspect()),
)
Colorbar(fig[:, end + 1], hm)

Overwrite a layer

There is an interpolate! method that will perform interpolation on an already existing layer. This is useful if you want to rapidly bring data from a layer to something compatible with another layer.

We can also pass a custom PROJ string if we want. For example, we will use an orthoraphic projection, centered around a point within the layer:

ws2 = interpolate(layer; dest = "+proj=ortho +lon_0=1.0859 +lat_0=45.6429")

🗺️  A 104 × 123 layer with 8556 Float32 cells
   Spatial Reference System: +proj=ortho +lat_0=45.6429 +lon_0=1 ... _defs

We can plot the outcome to see that it was correctly applied:

Code for the figure

fig, ax, hm = heatmap(
    ws2;
    colormap = :navia,
    figure = (; size = (800, 400)),
    axis = (; aspect = DataAspect()),
)
Colorbar(fig[:, end + 1], hm)

Related documentation

SimpleSDMLayers.interpolate Function

interpolate(layer::SDMLayer; dest="+proj=natearth2", newsize=nothing)

Returns an interpolated version of the later under the new destination CRS (natearth2 by default), and with optionally a new size of newsize.

source

interpolate(layer::SDMLayer, destination::SDMLayer)

Interpolates a layer target so that it uses the same grid, crs, etc as destination.

source

SimpleSDMLayers.interpolate! Function

interpolate!(destination::SDMLayer, source::SDMLayer)

Interpolates the data from the source layer into the destination layer. This is useful when the two layers have different resolutions, or different projections.

source

SimpleSDMLayers.projection Function

projection(layer::SDMLayer)

Returns the ArchGDAL SRS for a given layer.

source

SpeciesDistributionToolkit.reproject Function

reproject(pol::SimpleSDMPolygons.FeatureCollection, dest)

Converts the points in a polygon to a different CRS, defined by its WKT.

source

reproject(occ::AbstractOccurrenceCollection, dest)

Converts the points in an abstract occurrence collection to a different CRS, defined by its WKT. Note that this only returns the coordinates, and not the rest of the occurrence information.

source

reproject(sdm::AbstractSDM, dest)

Converts the georeferences instances in an SDM to a different CRS, defined by its WKT. Note that this only returns the coordinates, and not the rest of the model.

source