Skip to content

Commit

Permalink
Maybe breaking? use StepRangeLen instead of LinRange (#825)
Browse files Browse the repository at this point in the history
* use StepRangeLen instead of LinRange

* xsize

* length

* use range in resample tests

* fix resample range, segfault bug too
  • Loading branch information
rafaqz authored Dec 8, 2024
1 parent 32c2004 commit 11c6a6c
Show file tree
Hide file tree
Showing 4 changed files with 38 additions and 31 deletions.
4 changes: 2 additions & 2 deletions ext/RastersArchGDALExt/gdal_source.jl
Original file line number Diff line number Diff line change
Expand Up @@ -171,8 +171,8 @@ function RA._dims(raster::AG.RasterDataset, crs=nokw, mappedcrs=nokw)
xorder = xstep > 0 ? ForwardOrdered() : ReverseOrdered()
yorder = ystep > 0 ? ForwardOrdered() : ReverseOrdered()
# Create lookup index. LinRange is easiest always the right size after fp error
xindex = LinRange(xmin, xmax, xsize)
yindex = LinRange(ymax, ymin, ysize)
xindex = range(; start=xmin, stop=xmax, length=xsize)
yindex = range(; start=ymax, stop=ymin, length=ysize)

# Define `Projected` lookups fo X and Y dimensions
xlookup = Projected(xindex;
Expand Down
4 changes: 2 additions & 2 deletions src/sources/grd.jl
Original file line number Diff line number Diff line change
Expand Up @@ -74,8 +74,8 @@ function _dims(A::RasterDiskArray{GRDsource}, crs=nokw, mappedcrs=nokw)
# Not fully implemented yet
xy_metadata = _metadatadict(GRDsource())

xindex = LinRange(xbounds[1], xbounds[2] - xspan, xsize)
yindex = LinRange(ybounds[2] + yspan, ybounds[1], ysize)
xindex = range(; start=xbounds[1], stop=xbounds[2] - xspan, length=xsize)
yindex = range(; start=ybounds[2] + yspan, stop=ybounds[1], length=ysize)

xlookup = Projected(xindex;
order=GRD_X_ORDER,
Expand Down
59 changes: 33 additions & 26 deletions test/resample.jl
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
using Rasters, ArchGDAL, GeoInterface, Extents
using Rasters, ArchGDAL, GeoInterface
using Extents
using Test
using Rasters.Lookups
import DimensionalData: @dim, YDim
Expand Down Expand Up @@ -135,31 +136,6 @@ include(joinpath(dirname(pathof(Rasters)), "../test/test_utils.jl"))
@test eltype(r1) == UInt8
end

@testset "resample to the same size or resolution leaves raster unchanged" begin
res = 2
ys = LinRange((90 - res / 2):-res:(-90 + res / 2))
xs = LinRange((-180 + res / 2):res:(180 - res / 2))

point_dims = Y(ys), X(xs)
interval_dims = Y(ys; sampling=Intervals(Center())), X(xs; sampling=Intervals(Center()))
rev_y_point_dims = Y(reverse(ys)), X(xs)
rev_y_interval_dims = reverse(interval_dims[1]), interval_dims[2]
rev_x_point_dims = Y(ys), X(reverse(xs))
rev_x_interval_dims = interval_dims[1], reverse(interval_dims[2])
test_dims = (point_dims, interval_dims, rev_x_point_dims, rev_x_interval_dims, rev_y_point_dims, rev_y_interval_dims)
ds_fwd = point_dims; f = identity
ds_fwd = point_dims; f = reverse

for ds_fwd in test_dims, f in (identity, reverse)
ds = f(ds_fwd)
raster = Raster(rand(ds), crs=EPSG(4326), missingval=NaN)
resampled_res = resample(raster; res)
resampled_size = resample(raster; size=size(raster), method=:near)
@test resampled_res == resampled_size == raster
@test dims(resampled_res) == dims(resampled_size) == dims(raster)
end
end

@testset "dimensions matcha after resampling with only `to`" begin
# some weird dimensions
to = (
Expand All @@ -176,4 +152,35 @@ include(joinpath(dirname(pathof(Rasters)), "../test/test_utils.jl"))
@test dims(resampled_3D, (1,2)) == to
@test dims(resampled_3D, Z) == Z(1:2)
end

@testset "resample to the same size or resolution leaves raster unchanged" begin
res = 2
ys = range(; start=(90 - res / 2), step=-res, stop=(-90 + res / 2))
xs = range(; start=(-180 + res / 2), step=res, stop=(180 - res / 2))

point_dims = Y(ys), X(xs)
interval_dims = Y(ys; sampling=Intervals(Center())), X(xs; sampling=Intervals(Center()))
rev_y_point_dims = Y(reverse(ys)), X(xs)
rev_y_interval_dims = reverse(interval_dims[1]), interval_dims[2]
rev_x_point_dims = Y(ys), X(reverse(xs))
rev_x_interval_dims = interval_dims[1], reverse(interval_dims[2])
test_dims = (point_dims, interval_dims, rev_x_point_dims, rev_x_interval_dims, rev_y_point_dims, rev_y_interval_dims)
ds_fwd = point_dims; f = identity
ds_fwd = point_dims; f = reverse

ds = ds_fwd
raster = Raster(rand(ds), crs=EPSG(4326), missingval=NaN)
resampled_res = resample(raster; res)
resampled_size = resample(raster; size=size(raster), method=:near)
@test resampled_res == resampled_size == raster
@test dims(resampled_res) == dims(resampled_size) == dims(raster)

ds = reverse(ds_fwd)
raster = Raster(rand(ds), crs=EPSG(4326), missingval=NaN)
resampled_res = resample(raster; res)
resampled_size = resample(raster; size=size(raster), method=:near)
@test resampled_res == resampled_size == raster
@test dims(resampled_res) == dims(resampled_size) == dims(raster)
end

end
2 changes: 1 addition & 1 deletion test/sources/gdal.jl
Original file line number Diff line number Diff line change
Expand Up @@ -551,7 +551,7 @@ gdalpath = maybedownload(url)
@test order(dims(rast)) == (ForwardOrdered(), ForwardOrdered())
@test span(rast) == (Regular(1.0), Regular(1.0))
@test sampling(rast) == (Intervals(Start()), Intervals(Start()))
@test index(rast) == (LinRange(0.0, 239.0, 240), LinRange(0.0, 179.0, 180))
@test index(rast) == (range(; start=0.0, stop=239.0, length=240), range(start=0.0, stop=179.0, length=180))
end

end
Expand Down

0 comments on commit 11c6a6c

Please sign in to comment.