rafaqz / DimensionalData.jl

Named dimensions and indexing for julia arrays and other data
https://rafaqz.github.io/DimensionalData.jl/stable/
MIT License
269 stars 37 forks source link

Indexing with a point extent returns zero dimensional array #578

Closed felixcremer closed 6 months ago

felixcremer commented 8 months ago

When I try to subset a Raster or YAXArray with an Extent which I got from a point and therefore the elements are the same I get a zero dimensional array back. I would have expected to get the pixel in which the point lays.

This is a reduced example using only DimArray. We could discuss whether this should be only possible to do for Intervals In my real life examples the X and Y dimensions are Intervals and there I would expect that I get the Interval that the extent overlays. It would also nice to be able do arr[Near(Extent(X=(1,1), Y=(2.1,2.1))] to specify that I want the nearest point to my extent.

julia> arr = DimArray(reshape(1:100,10,10), (X(1:10), Y(1:10)))
10×10 DimArray{Int64,2} with dimensions: 
  X Sampled{Int64} 1:10 ForwardOrdered Regular Points,
  Y Sampled{Int64} 1:10 ForwardOrdered Regular Points
      1   2   3   4   5   6   7   8   9   10
  1   1  11  21  31  41  51  61  71  81   91
  2   2  12  22  32  42  52  62  72  82   92
  ⋮                   ⋮                    ⋮
  8   8  18  28  38  48  58  68  78  88   98
  9   9  19  29  39  49  59  69  79  89   99
 10  10  20  30  40  50  60  70  80  90  100

julia> arr[Extent(X=(1,1), Y=(2.1,2.1))]
1×0 DimArray{Int64,2} with dimensions: 
  X Sampled{Int64} 1:1 ForwardOrdered Regular Points,
  Y Sampled{Int64} 3:2 ForwardOrdered Regular Points

 1
julia> arr[Extent(X=(1,1), Y=(2,2))] # This works, because it is the exact values in the dimensions.
1×1 DimArray{Int64,2} with dimensions: 
  X Sampled{Int64} 1:1 ForwardOrdered Regular Points,
  Y Sampled{Int64} 2:2 ForwardOrdered Regular Points
     2
 1  11
rafaqz commented 8 months ago

I think you can wrap it with Touches(extent) to get the one pixel? Selecting with an extent with no area would only return an exact point from a Points raster.

felixcremer commented 8 months ago

Wrapping it into Touches works for Rasters but not for YAXArray. With a YAXArray DimensionalData gives a range that is of by one and thereby I get an empty range. I suspect that this is coming from the fact, that The dimensions in the YAXArray are all point axes. the upper one is the dimensions of the YAXArray and below are the dimensions of the Raster:

julia> dims(rxy)
X Sampled{Float64} -28.4029998779297:0.0550011811599821:18.18300056457514 ForwardOrdered Regular Points,
Y Sampled{Float64} -23.4029998779297:0.0550012159753107:21.863000869751005 ForwardOrdered Regular Points,
Ti Sampled{DateTime} DateTime[1996-01-16T12:00:00, …, 1996-12-16T12:00:00] ForwardOrdered Irregular Points

julia> dims(rasxy)
X Sampled{Float64} Float64[-28.4029998779297, -28.347998696769718, …, 18.12799938341516, 18.18300056457514] ForwardOrdered Explicit Intervals{Center},
Y Sampled{Float64} Float64[-23.4029998779297, -23.347998661954392, …, 21.807999653775695, 21.863000869751005] ForwardOrdered Explicit Intervals{Center},
Ti Sampled{DateTime} DateTime[1996-01-16T12:00:00, …, 1996-12-16T12:00:00] ForwardOrdered Explicit Intervals{Center}

Does this mean, that we would have to ensure that the dimensions are Intervals in YAXArrays?

rafaqz commented 8 months ago

Im not even sure what it does for points or why its different! It could be a bug, I probably assumed Touches was only for Intervals.

But it does sound like what you have is interval pixels if you want a point to fall somewhere inside a pixel.

You can detect if a netcdf or tif is intervals, and most of the time it usually is.

rafaqz commented 6 months ago

Think I will close this, its really just YAXArrays using Points when it has Intervals