... only read part of a layer?

using SpeciesDistributionToolkit
using CairoMakie

To illustrate the process, we will load the demo dataset:

t = SimpleSDMLayers.__demodata()
f = tempname() * ".tiff"
SimpleSDMLayers.save(f, t)

"/tmp/jl_9G3YJMN6D4.tiff"

We can check the look of this layer as a reference:

Code for the figure

heatmap(t; colormap = :navia)

We now define a bounding box:

bbox = (left = -76.0, right = -73.0, bottom = 48.0, top = 50.0)

(left = -76.0, right = -73.0, bottom = 48.0, top = 50.0)

WGS84 is used for the bounding box

We are assuming that there is not information about the CRS of the layer we are working on until it is read, and for this reason the bounding box limits are given in lat/lon.

We can now pass all of these arguments to the SDMLayer function:

k = SDMLayer(f; bandnumber = 1, bbox...)

SDM Layer with 57354 UInt16 cells
	Proj string: +proj=aea +lat_0=44 +lon_0=-68.5 +lat_1=60 +lat_2=46 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs
	Grid size: (242, 237)

Note that this layer is indeed cropped, but has retained its CRS:

Code for the figure

heatmap(k; colormap = :navia)

We can also plot this layer projected under WGS84, and overlay the boundingbox:

poly = [
    (bbox.left, bbox.bottom),
    (bbox.right, bbox.bottom),
    (bbox.right, bbox.top),
    (bbox.left, bbox.top),
    (bbox.left, bbox.bottom),
]

5-element Vector{Tuple{Float64, Float64}}:
 (-76.0, 48.0)
 (-73.0, 48.0)
 (-73.0, 50.0)
 (-76.0, 50.0)
 (-76.0, 48.0)

Code for the figure

heatmap(
    interpolate(k; dest = "+proj=longlat +datum=WGS84 +no_defs +type=crs", newsize=(500, 500));
    colormap = :navia,
)
lines!(poly; color = :black, linewidth = 2, linestyle = :dash)

Note that the definition of the bounding box is that all of the content in the bounding box is included in the final layer, even if it results in some values outside the bounding box being included as well.