Rasters.jl
Rasters.jl copied to clipboard
Negative indexing does not reverse Dims
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.
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.
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.
By the way the easy fix is:
A = set(A, Y => DimentionalData.ReverseOrdered())