From 670efcd39545d6131f6666e373c96510f5d767d6 Mon Sep 17 00:00:00 2001
From: tpoisot <timothee.poisot@umontreal.ca>
Date: Thu, 19 Sep 2024 20:57:03 -0400
Subject: [PATCH] bug(layers): The bounding box should be converted from WGS84
 to the layer CRS #270

---
 SimpleSDMLayers/Project.toml         |  2 +-
 SimpleSDMLayers/src/io/geotiff.jl    |  6 ++++++
 SimpleSDMLayers/src/io/read_write.jl | 16 ++++++++++++++++
 3 files changed, 23 insertions(+), 1 deletion(-)

diff --git a/SimpleSDMLayers/Project.toml b/SimpleSDMLayers/Project.toml
index 9f2064917..58f303ff4 100644
--- a/SimpleSDMLayers/Project.toml
+++ b/SimpleSDMLayers/Project.toml
@@ -1,7 +1,7 @@
 name = "SimpleSDMLayers"
 uuid = "2c645270-77db-11e9-22c3-0f302a89c64c"
 authors = ["Timothée Poisot <timothee.poisot@umontreal.ca>", "Gabriel Dansereau <gabriel.dansereau@umontreal.ca>"]
-version = "1.0.3"
+version = "1.0.4"
 
 [deps]
 ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"
diff --git a/SimpleSDMLayers/src/io/geotiff.jl b/SimpleSDMLayers/src/io/geotiff.jl
index 3929be770..e1c55ab73 100644
--- a/SimpleSDMLayers/src/io/geotiff.jl
+++ b/SimpleSDMLayers/src/io/geotiff.jl
@@ -71,6 +71,12 @@ function _read_geotiff(
         maxlon = minlon + size(band, 1) * transform[2]
         minlat = maxlat - abs(size(band, 2) * transform[6])
 
+        # We need to make sure the WGS84 coordinates for the boundingbox are included in the layer
+        prj = Proj.Transformation("+proj=longlat +datum=WGS84 +no_defs +type=crs", wkt; always_xy=true)
+        left, bottom = prj(left, bottom)
+        right, top = prj(right, top)
+        
+        # And now we crop
         left = isnothing(left) ? minlon : max(left, minlon)
         right = isnothing(right) ? maxlon : min(right, maxlon)
         bottom = isnothing(bottom) ? minlat : max(bottom, minlat)
diff --git a/SimpleSDMLayers/src/io/read_write.jl b/SimpleSDMLayers/src/io/read_write.jl
index 8fd42acff..d187160b6 100644
--- a/SimpleSDMLayers/src/io/read_write.jl
+++ b/SimpleSDMLayers/src/io/read_write.jl
@@ -37,4 +37,20 @@ end
     SimpleSDMLayers.save(f, t)
     k = SDMLayer(f)
     @test SimpleSDMLayers._layers_are_compatible(t, k)
+end
+
+@testitem "We can save a layer and read with the correct bbox" begin
+    t = SimpleSDMLayers.__demodata(; reduced=true)
+    f = tempname()*".tiff"
+    f = "test.tiff"
+    SimpleSDMLayers.save(f, t)
+    # WGS84 smaller bounding box
+    bbox = (left=-79., right=-75., bottom=47., top=49.)
+    k = SDMLayer(f; bandnumber=1, bbox...)
+    @test all(size(k) .< size(t))
+    @test k.crs == t.crs
+    @test k.x[1] > t.x[1]
+    @test k.x[2] < t.x[2]
+    @test k.y[1] > t.y[1]
+    @test k.y[2] < t.y[2]
 end
\ No newline at end of file