From d56f46b311f0707b136313ace78dc84f93e99c1f Mon Sep 17 00:00:00 2001
From: Gabriel Dansereau <gabriel.dansereau@umontreal.ca>
Date: Fri, 4 Jun 2021 12:04:01 -0400
Subject: [PATCH] =?UTF-8?q?=F0=9F=90=9B=20fix=20geotiff=20writing=20with?=
 =?UTF-8?q?=20custom=20nodata=20field,=20closes=20#108?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 src/datasets/geotiff.jl | 10 +++++-----
 test/dataread.jl        | 20 +++++++++++++++-----
 test/subsetting.jl      |  2 +-
 3 files changed, 21 insertions(+), 11 deletions(-)

diff --git a/src/datasets/geotiff.jl b/src/datasets/geotiff.jl
index e84f30f0..1694ed72 100644
--- a/src/datasets/geotiff.jl
+++ b/src/datasets/geotiff.jl
@@ -102,7 +102,7 @@ Write a single `layer` to a `file`, where the `nodata` field is set to an
 arbitrary value.
 """
 function geotiff(file::AbstractString, layer::SimpleSDMPredictor{T}; nodata::T=convert(T, -9999)) where {T <: Number}
-    array_t = _prepare_layer_for_burnin(layer)
+    array_t = _prepare_layer_for_burnin(layer, nodata)
     width, height = size(array_t)
 
     # Geotransform
@@ -138,9 +138,9 @@ function geotiff(file::AbstractString, layer::SimpleSDMPredictor{T}; nodata::T=c
     return file
 end
 
-function _prepare_layer_for_burnin(layer::T) where {T <: SimpleSDMLayer}
+function _prepare_layer_for_burnin(layer::SimpleSDMPredictor{T}, nodata::T) where {T <: Number}
     @assert eltype(layer) <: Number
-    array = replace(layer.grid, nothing => NaN)
+    array = replace(layer.grid, nothing => nodata)
     array = convert(Matrix{eltype(layer)}, array)
     dtype = eltype(array)
     array_t = reverse(permutedims(array, [2, 1]); dims=2)
@@ -156,7 +156,7 @@ Stores a series of `layers` in a `file`, where every layer in a band. See
 function geotiff(file::AbstractString, layers::Vector{SimpleSDMPredictor{T}}; nodata::T=convert(T, -9999)) where {T <: Number}
     bands = 1:length(layers)
     _layers_are_compatible(layers)
-    width, height = size(_prepare_layer_for_burnin(layers[1]))
+    width, height = size(_prepare_layer_for_burnin(layers[1], nodata))
 
     # Geotransform
     gt = zeros(Float64, 6)
@@ -179,7 +179,7 @@ function geotiff(file::AbstractString, layers::Vector{SimpleSDMPredictor{T}}; no
             band = ArchGDAL.getband(dataset, i)
             
             # Write data to band
-            ArchGDAL.write!(band, _prepare_layer_for_burnin(layers[i]))
+            ArchGDAL.write!(band, _prepare_layer_for_burnin(layers[i], nodata))
 
             # Write nodata and projection info
             ArchGDAL.setnodatavalue!(band, nodata)
diff --git a/test/dataread.jl b/test/dataread.jl
index df4df3c7..48ef2692 100644
--- a/test/dataread.jl
+++ b/test/dataread.jl
@@ -3,11 +3,21 @@ using SimpleSDMLayers
 using Test
 
 l = SimpleSDMPredictor(WorldClim, BioClim, 1)
-f = tempname()
-geotiff(f, l)
-mp = geotiff(SimpleSDMResponse, f)
 
-@test typeof(mp) <: SimpleSDMResponse
-@test size(mp) == size(l)
+f1 = tempname()
+geotiff(f1, l)
+mp1 = geotiff(SimpleSDMResponse, f1)
+
+@test typeof(mp1) <: SimpleSDMResponse
+@test size(mp1) == size(l)
+@test mp1 == l
+
+f2 = tempname()
+geotiff(f2, l; nodata=-3.4f38)
+mp2 = geotiff(SimpleSDMPredictor, f2)
+
+@test typeof(mp2) <: SimpleSDMResponse
+@test size(mp2) == size(l)
+@test mp2 == l
 
 end
diff --git a/test/subsetting.jl b/test/subsetting.jl
index bcbed4b2..7f16bb85 100644
--- a/test/subsetting.jl
+++ b/test/subsetting.jl
@@ -10,7 +10,7 @@ l1 = temp[coords]
 l2 = SimpleSDMPredictor(WorldClim, BioClim, 1; coords...)
 tempfile = tempname()
 geotiff(tempfile, l2)
-l3 = replace(geotiff(SimpleSDMPredictor, tempfile), NaN => nothing)
+l3 = geotiff(SimpleSDMPredictor, tempfile)
 
 @test size(l1) == size(l2)
 @test size(l1) == size(l3)