gabrieldansereau commented 3 years ago

What the pull request does

:boom: Breaking changes! :bug: But it's actually just solving bugs and unexpected behaviours, some of which have been around since v0.4.3. The calls should work as before but will return different results in important functions.

Fixes #99, fixes #100, fixes #102, fixes #103, fixes #104.

Checklist

Using:

layer = SimpleSDMPredictor(WorldClim, BioClim, 1)
coords = (left = -145.0, right = -50.0, bottom = 20.0, top = 75.0)

This PR needs to fix:

[x] getindex(layer; coords) should return a layer with dimensions 330 x 570
[x] getindex(layer; coords) should return a layer whose bounding coordinates are exactly (20.083333333333332, 74.91666666666667) for latitude and (-144.91666666666666, -50.083333333333336) for longitude (same as 10 arc-min WorldClim raster cropped in QGIS/GDAL)
[x] getindex(layer; coords) should return a layer whose stride is exactly (0.08333333333333333, 0.08333333333333333)
[x] geotiff(layer; coords...) should return the same thing
[x] geotiff("test.tif", layer) and geotiff(SimpleSDMPredictor, "test.tif") should allow to write and re-read the same layer
[x] All these changes should have appropriate tests

See my before/after comment for the details.

Type of change

Please indicate the relevant option(s)

[x] :bug: Bug fix (non-breaking change which fixes an issue)
[ ] :sparkle: New feature (non-breaking change which adds functionality)
[x] :boom: Breaking change (fix or feature that would cause existing functionality to not work as expected)
[x] :book: This change requires a documentation update

Checklist

[x] The changes are documented
- [x] The docstrings of the different functions describe the arguments, purpose, and behavior
- [x] There are examples in the documentation website
[x] The changes are tested
[x] The changes do not modify the behavior of the previously existing functions
- If they do, please explain why and how in the introduction paragraph
[x] For new contributors - my name and information are added to .zenodo.json
[x] The Project.toml field version has been updated
- Change the last number for a v0.0.x series release, a bug fix, or a performance improvement
- Change the middle number for additional features that do not break the current test suite (unless you fix a bug in the test suite)
- Change the first number for changes that break the current test suite

codecov-commenter commented 3 years ago

Codecov Report

Merging #101 (9423e93) into master (f27ea63) will increase coverage by 1.52%. The diff coverage is 96.93%.

@@            Coverage Diff             @@
##           master     #101      +/-   ##
==========================================
+ Coverage   80.02%   81.54%   +1.52%     
==========================================
  Files          24       24              
  Lines         756      775      +19     
==========================================
+ Hits          605      632      +27     
+ Misses        151      143       -8

Flag	Coverage Δ
unittests	`81.54% <96.93%> (+1.52%)`	:arrow_up:

Flags with carried forward coverage won't be shown. Click here to find out more.

Impacted Files	Coverage Δ
src/lib/overloads.jl	`88.66% <96.05%> (+1.74%)`	:arrow_up:
src/datasets/geotiff.jl	`73.68% <100.00%> (+9.55%)`	:arrow_up:
src/lib/basics.jl	`100.00% <100.00%> (ø)`
src/lib/types.jl	`100.00% <100.00%> (ø)`
src/SimpleSDMLayers.jl	`80.00% <0.00%> (-12.31%)`	:arrow_down:

Continue to review full report at Codecov.

Legend - Click here to learn more Δ = absolute <relative> (impact), ø = not affected, ? = missing data Powered by Codecov. Last update f27ea63...9423e93. Read the comment docs.

gabrieldansereau commented 3 years ago

Here's a before & after of the changes. Using this:

temp = SimpleSDMPredictor(WorldClim, BioClim, 1);
coords = (left = -145.0, right = -50.0, bottom = 20.0, top = 75.0);

# getindex
l1 = temp[coords];
size(l1)
(l1.bottom, l1.top)
(l1.left, l1.right)
stride(l1)

# geotiff reading
l2 = SimpleSDMPredictor(WorldClim, BioClim, 1; coords...);
l2.grid == l1.grid
size(l2)
(l2.bottom, l2.top)
(l2.left, l2.right)
stride(l2)

# geotiff writing
l2_backup = copy(l2);
tempfile = tempname();
geotiff(tempfile, l2);
l3 = replace(geotiff(SimpleSDMPredictor, tempfile), NaN => nothing);
l2_backup.grid == l2.grid
l2_backup.grid == l3.grid

# Inverse bounds
SimpleSDMPredictor(temp.grid, 180.0, -180.0, 90.0, -90.0)

# hcat/vcat
l4 = SimpleSDMPredictor(WorldClim, BioClim, 1; left=0.0, right=10.0, bottom=0.0, top=10.0);
l5 = SimpleSDMPredictor(WorldClim, BioClim, 1; left=0.0, right=10.0, bottom=10.0, top=20.0);
l6 = SimpleSDMPredictor(WorldClim, BioClim, 1; left=0.0, right=10.0, bottom=20.0, top=30.0);
vcat(l4, l5)
vcat(l4, l6)

Before:

# getindex
julia> l1 = temp[coords];
julia> size(l1)
(332, 571)
julia> (l1.bottom, l1.top)
(19.83333333333334, 75.16666666666666)
julia> (l1.left, l1.right)
(-145.16666666666669, -49.99999999999998)
julia> stride(l1)
(0.08333333333333337, 0.0833333333333333)

# geotiff reading
julia> l2 = SimpleSDMPredictor(WorldClim, BioClim, 1; coords...);
julia> l2.grid == l1.grid
false
julia> size(l2)
(331, 571)
julia> (l2.bottom, l2.top)
(19.833333333333336, 75.0)
julia> (l2.left, l2.right)
(-145.16666666666669, -50.0)
julia> stride(l2)
(0.08333333333333336, 0.08333333333333333)

# geotiff writing
julia> l2_backup = copy(l2);
julia> tempfile = tempname();
julia> geotiff(tempfile, l2);
julia> l3 = replace(geotiff(SimpleSDMPredictor, tempfile), NaN => nothing);
julia> l2_backup.grid == l2.grid
false
julia> l2_backup.grid == l3.grid
false

# Inverse bounds
julia> SimpleSDMPredictor(temp.grid, 180.0, -180.0, 90.0, -90.0)
SDM predictor → 1080×2160 grid with 808053 Float32-valued cells
  Latitudes     (-89.91666666666667, 89.91666666666667)
  Longitudes    (-179.91666666666666, 179.91666666666666)

# hcat/vcat
julia> l4 = SimpleSDMPredictor(WorldClim, BioClim, 1; left=0.0, right=10.0, bottom=0.0, top=10.0);
julia> l5 = SimpleSDMPredictor(WorldClim, BioClim, 1; left=0.0, right=10.0, bottom=10.0, top=20.0);
julia> l6 = SimpleSDMPredictor(WorldClim, BioClim, 1; left=0.0, right=10.0, bottom=20.0, top=30.0);
julia> vcat(l4, l5)
SDM predictor → 122×61 grid with 5573 Float32-valued cells
  Latitudes     (9.834016393442623, 9.99931693989071)
  Longitudes    (-0.08333333333333333, 9.916666666666666)
julia> vcat(l4, l6)
SDM predictor → 122×61 grid with 5573 Float32-valued cells
  Latitudes     (10.040300546448087, 19.79303278688525)
  Longitudes    (-0.08333333333333333, 9.916666666666666)

After:

# getindex
julia> l1 = temp[coords];
julia> size(l1)
(330, 570)
julia> (l1.bottom, l1.top)
(20.0, 75.0)
julia> (l1.left, l1.right)
(-145.0, -50.0)
julia> stride(l1)
(0.08333333333333333, 0.08333333333333333)

# geotiff reading
julia> l2 = SimpleSDMPredictor(WorldClim, BioClim, 1; coords...);
julia> l2.grid == l1.grid
true
julia> size(l2)
(330, 570)
julia> (l2.bottom, l2.top)
(20.0, 75.0)
julia> (l2.left, l2.right)
(-145.0, -50.0)
julia> stride(l2)
(0.08333333333333333, 0.08333333333333333)

# geotiff writing
julia>  l2_backup = copy(l2);
julia> tempfile = tempname();
julia> geotiff(tempfile, l2);
julia> l3 = replace(geotiff(SimpleSDMPredictor, tempfile), NaN => nothing);
julia> l2_backup.grid == l2.grid
true
julia> l2_backup.grid == l3.grid
true

# Inverse bounds
julia> SimpleSDMPredictor(temp.grid, 180.0, -180.0, 90.0, -90.0)
ERROR: LoadError: ArgumentError: The right bounding coordinate must be greater than the right one

# hcat/vcat
julia> l4 = SimpleSDMPredictor(WorldClim, BioClim, 1; left=0.0, right=10.0, bottom=0.0, top=10.0);
julia> l5 = SimpleSDMPredictor(WorldClim, BioClim, 1; left=0.0, right=10.0, bottom=10.0, top=20.0);
julia> l6 = SimpleSDMPredictor(WorldClim, BioClim, 1; left=0.0, right=10.0, bottom=20.0, top=30.0);
julia> vcat(l4, l5)
SDM predictor → 120×60 grid with 5419 Float32-valued cells
  Latitudes     (0.08333333333333333, 19.916666666666668)
  Longitudes    (0.08333333333333333, 9.916666666666666)
julia> vcat(l4, l6)
ERROR: LoadError: ArgumentError: The two layers passed to vcat must have contiguous bottom and top coordinates

gabrieldansereau commented 3 years ago

@tpoisot This fixes the bugs I've had for a while.

A couple of things you might have a different opinion on:

I needed to switch back to range instead of LinRange, otherwise I couldn't find a proper fix for one of my problems. As I understand, since we're doing operations such as longitudes(layer)+stride(layer,1), range is preferable as it will avoid approximation errors. I didn't seem to add that many allocations and made many things simpler.
My fix for getindex() and geotiff(): I added a side=:left keyword argument in _match_longitude() et al. which solves the dimensions problem.
I left the default as side=:none so it doesn't break the behaviour in other functions. However, this might create an ambiguous situation when getting a coordinate that is exactly between two cells. layer[lon1, lat1] might return a different value than layer[lon1:lon2, lat1:lat2][1] for instance.
Layers with inversed coordinates aren't possible anymore (#103)
hcat and vcat returned super weird results (#102 ). I think I fixed it, unless that behaviour was on purpose.

tpoisot commented 3 years ago

feel free to merge and tag the new release - good job

PoisotLab / SimpleSDMLayers.jl

Fix layer dimensions bugs #101

Codecov Report