Rasters.jl icon indicating copy to clipboard operation
Rasters.jl copied to clipboard

Negative indexing does not reverse Dims

Open mauro3 opened this issue 1 year ago • 3 comments

using Rasters, Plots
url = "https://www.unidata.ucar.edu/software/netcdf/examples/tos_O1_2001-2002.nc";
filename = download(url, "tos_O1_2001-2002.nc");
A = Raster(filename)

now runing

julia> dims(A[:,end:-1:1,:])
X Mapped Float64[1.0, 3.0, …, 357.0, 359.0] ForwardOrdered Explicit Intervals crs: EPSG mappedcrs: EPSG,
Y Mapped Float64[89.5, 88.5, …, -78.5, -79.5] ForwardOrdered Explicit Intervals crs: EPSG mappedcrs: EPSG,
Ti Sampled CFTime.DateTime360Day[CFTime.DateTime360Day(2001-01-16T00:00:00), …, CFTime.DateTime360Day(2002-12-16T00:00:00)] ForwardOrdered Explicit Intervals

whereas

julia> dims(reverse(A, dims=2))
X Mapped Float64[1.0, 3.0, …, 357.0, 359.0] ForwardOrdered Explicit Intervals crs: EPSG mappedcrs: EPSG,
Y Mapped Float64[89.5, 88.5, …, -78.5, -79.5] ReverseOrdered Explicit Intervals crs: EPSG mappedcrs: EPSG,
Ti Sampled CFTime.DateTime360Day[CFTime.DateTime360Day(2001-01-16T00:00:00), …, CFTime.DateTime360Day(2002-12-16T00:00:00)] ForwardOrdered Explicit Intervals

i.e. once Y is ForwardOrdered once ReverseOrdered. Shouldn't they be the same? Using negative indexing leads to issues, for instance, with plotting.

mauro3 avatar Sep 08 '22 09:09 mauro3

Totally, there are some problems with unorthodox indexing like this in DimensionalData.jl, because it relies on searchsorted* so much.

The problem is there is no entirely type-stable way to handle a reverse index and have also fast selector lookups. searchsorted is fastest when the sort direction is known at compile time. So we store ForwardOrdered or ReverseOrdered in the type. But that means changing the order by indexing in reverse has undefined behavior. So using reverse is the preferred method.

We probably just need to check all StepRange and similar ranges and make indexing with them type unstable for all ordered lookups.

rafaqz avatar Sep 08 '22 10:09 rafaqz

Annoyingly, updating a LookupArray when indexing with a steprange goes from 17ns to 250ns with this instability. It's nothing when you are working with large arrays but annoying if you're using very small arrays or just taking a view. It also means taking a view of a Raster with a StepRange wont work in a GPU kernel or anywhere else you need type stability. Although I guess a UnitRange is more often used for that.

Maybe it should just throw an error instead and tell you to use reverse. Its the only way I can see that keeps type stability and doesn't let the weird behaviour in.

rafaqz avatar Sep 08 '22 20:09 rafaqz

By the way the easy fix is:

A = set(A, Y => DimentionalData.ReverseOrdered())

rafaqz avatar Sep 13 '22 22:09 rafaqz