precise step (#656)

* calculate step more precisely and maybe round it

* do not do try to improve accuracy for integers

* Update src/sources/commondatamodel.jl

* rounding atol is based on index eltype

* add tests

---------

Co-authored-by: Rafael Schouten <[email protected]>

src/sources/commondatamodel.jl

@@ -385,7 +385,15 @@ end
 function _cdmspan(index, order)
     # Handle a length 1 index
     length(index) == 1 && return Regular(zero(eltype(index))), Points()
-    step = index[2] - index[1]
+
+    step = if eltype(index) <: AbstractFloat
+        # Calculate step, avoiding as many floating point errors as possible
+        st = Base.step(Base.range(Float64(first(index)), Float64(last(index)); length = length(index)))
+        st_rd = round(st, digits = Base.floor(Int,-log10(eps(eltype(index))))) # round to nearest digit within machine epsilon
+        isapprox(st_rd, st; atol = eps(eltype(index))) ? st_rd : st # keep the rounded number if it is very close to the original
+    else
+        index[2] - index[1]
+    end
     for i in 2:length(index)-1
         # If any step sizes don't match, its Irregular
         if !(index[i+1] - index[i] ≈ step)

test/sources/commondatamodel.jl

@@ -0,0 +1,24 @@
+using Rasters, NCDatasets, Test
+import Rasters: ForwardOrdered, ReverseOrdered, Regular
+@testset "step" begin
+    # test if regular indices are correctly rounded
+    f32_indices = range(0.075f0, 10.075f0; step = 0.05f0) |> collect
+    @test Rasters._cdmspan(f32_indices, ForwardOrdered())[1] === Regular(0.05)
+
+    f32_indices_rev = range(10.075f0, 0.075f0; step = -0.05f0) |> collect
+    @test Rasters._cdmspan(f32_indices_rev, ReverseOrdered())[1] === Regular(-0.05)
+
+    # test if regular indices are not rounded when they should not
+    indices_one_third = range(0, 10; length = 31) |> collect
+    @test Rasters._cdmspan(indices_one_third, ForwardOrdered())[1] === Regular(1/3)
+
+    # test when reading a file
+    ras = Raster(rand(X(f32_indices), Y(indices_one_third)))
+    tempfile = tempname() * ".nc"
+    write(tempfile, ras)
+    ras_read = Raster(tempfile)
+    steps = step.(dims(ras_read))
+    @test steps[1] == 0.05
+    @test steps[2] == 1/3
+
+end