bug(layers): The bounding box should be converted from WGS84 to the l…

…ayer CRS #270

SimpleSDMLayers/Project.toml

@@ -1,7 +1,7 @@
 name = "SimpleSDMLayers"
 uuid = "2c645270-77db-11e9-22c3-0f302a89c64c"
 authors = ["Timothée Poisot <[email protected]>", "Gabriel Dansereau <[email protected]>"]
-version = "1.0.3"
+version = "1.0.4"
 
 [deps]
 ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"

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)

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