WIP: geometry lookup #854

asinghvi17 · 2025-01-08T11:29:22Z

TODOs:

felixcremer · 2025-04-22T09:10:51Z

How far is this? Is there anything I could help to push this through?

asinghvi17 · 2025-04-22T11:00:23Z

The big things we need are:

Zonal using geometrylookup to return a vector data cube
Zonal on VDC to give a view of all geoms that intersect the input geoms
resample from one geometrylookup to another

All of these are easy enough. I can push those today. What i really need is IO support via CF standards. I can probably figure out how to output but no idea how to parse input. We also need the dimarrays to tables thing to parse GeoJSON to a vdc, or some way to convert a GeoJSON wide table to a dimarray.

rafaqz · 2025-04-22T12:33:19Z

Yeah really the big task we haven't started is having geometry writing for CF, maybe in CommonDataModel.jl

… wraps that

This way we can force skipmissing=false and manually do it, if we want e.g. vector data cubes.

asinghvi17 · 2025-04-22T15:41:51Z

src/methods/zonal.jl

@rafaqz would like to get your opinion on this, not sure how breaking this will be to end users, who may now get an empty array where previously the function they provided wasn't even called.

from Slack this is extremely cool and we should do this

asinghvi17 · 2025-04-22T15:46:33Z

I guess the place to do IO is:

_cdmlookup for read
_def_dim_var! for write
?

with style and ease

rafaqz · 2025-04-22T19:24:07Z

Maybe we can just implement it all here first if that's easier

asinghvi17 · 2025-04-22T19:45:28Z

Yeah that's what I'm trying to do with zonal for now, then I can just paste that file and tests to a new PR that does this for all zonal

so that we can check that the thing is explicitly false, this is quite useful for spatialslice because you can provide anything but false there, and it assumes true

we need an interface / spec for what rebuild must accept in DD!!

felixcremer · 2025-04-24T13:53:31Z

Near on a GeometryLookup should work but gives the following error:

Stacktrace:
  [1] _selnotfound(l::GeometryLookup{…}, selval::Float64)
    @ DimensionalData.Dimensions.Lookups ~/Documents/VectorDataCubes_Tutorial_EGU_2025/dev/DimensionalData/src/Lookups/selector.jl:233
  [2] _selnotfound_or_nothing(err::DimensionalData.Dimensions.Lookups._True, lookup::GeometryLookup{…}, selval::Float64)
    @ DimensionalData.Dimensions.Lookups ~/Documents/VectorDataCubes_Tutorial_EGU_2025/dev/DimensionalData/src/Lookups/selector.jl:231
  [3] at(::DimensionalData.Dimensions.Lookups.Unordered, ::DimensionalData.Dimensions.Lookups.NoSpan, lookup::GeometryLookup{…}, selval::Float64, atol::Nothing, rtol::Nothing; err::DimensionalData.Dimensions.Lookups._True)
    @ DimensionalData.Dimensions.Lookups ~/Documents/VectorDataCubes_Tutorial_EGU_2025/dev/DimensionalData/src/Lookups/selector.jl:216
  [4] at(::DimensionalData.Dimensions.Lookups.Unordered, ::DimensionalData.Dimensions.Lookups.NoSpan, lookup::GeometryLookup{…}, selval::Float64, atol::Nothing, rtol::Nothing)
    @ DimensionalData.Dimensions.Lookups ~/Documents/VectorDataCubes_Tutorial_EGU_2025/dev/DimensionalData/src/Lookups/selector.jl:213
  [5] at(lookup::GeometryLookup{…}, sel::At{…}; kw::@Kwargs{})
    @ DimensionalData.Dimensions.Lookups ~/Documents/VectorDataCubes_Tutorial_EGU_2025/dev/DimensionalData/src/Lookups/selector.jl:158
  [6] at(lookup::GeometryLookup{…}, sel::At{…})
    @ DimensionalData.Dimensions.Lookups ~/Documents/VectorDataCubes_Tutorial_EGU_2025/dev/DimensionalData/src/Lookups/selector.jl:157
  [7] near(order::DimensionalData.Dimensions.Lookups.Unordered, ::DimensionalData.Dimensions.Lookups.NoSampling, lookup::GeometryLookup{…}, sel::Near{…}; kw::@Kwargs{})
    @ DimensionalData.Dimensions.Lookups ~/Documents/VectorDataCubes_Tutorial_EGU_2025/dev/DimensionalData/src/Lookups/selector.jl:285
  [8] near(lookup::GeometryLookup{…}, sel::Near{…}; kw::@Kwargs{})
    @ DimensionalData.Dimensions.Lookups ~/Documents/VectorDataCubes_Tutorial_EGU_2025/dev/DimensionalData/src/Lookups/selector.jl:283
  [9] near(lookup::GeometryLookup{…}, sel::Near{…})
    @ DimensionalData.Dimensions.Lookups ~/Documents/VectorDataCubes_Tutorial_EGU_2025/dev/DimensionalData/src/Lookups/selector.jl:277
 [10] _selectindices
    @ ~/Documents/VectorDataCubes_Tutorial_EGU_2025/dev/DimensionalData/src/Lookups/selector.jl:1098 [inlined]
 [11] selectindices(l::GeometryLookup{…}, sel::Near{…}; kw::@Kwargs{})
    @ DimensionalData.Dimensions.Lookups ~/Documents/VectorDataCubes_Tutorial_EGU_2025/dev/DimensionalData/src/Lookups/selector.jl:1075
 [12] #_selecttuple#70
    @ ~/Documents/VectorDataCubes_Tutorial_EGU_2025/dev/DimensionalData/src/Lookups/selector.jl:370 [inlined]
 [13] _selecttuple
    @ ~/Documents/VectorDataCubes_Tutorial_EGU_2025/dev/DimensionalData/src/Lookups/selector.jl:368 [inlined]
 [14] selectindices
    @ ~/Documents/VectorDataCubes_Tutorial_EGU_2025/dev/DimensionalData/src/Lookups/selector.jl:1069 [inlined]
 [15] _dims2indices
    @ ~/Documents/VectorDataCubes_Tutorial_EGU_2025/dev/DimensionalData/src/Dimensions/indexing.jl:146 [inlined]
 [16] _dims2indices
    @ ~/Documents/VectorDataCubes_Tutorial_EGU_2025/dev/DimensionalData/src/Dimensions/indexing.jl:138 [inlined]
 [17] dims2indices
    @ ~/Documents/VectorDataCubes_Tutorial_EGU_2025/dev/DimensionalData/src/Dimensions/indexing.jl:31 [inlined]
 [18] macro expansion
    @ ~/Documents/VectorDataCubes_Tutorial_EGU_2025/dev/DimensionalData/src/Dimensions/indexing.jl:91 [inlined]
 [19] split_alignments(fa::typeof(DimensionalData.Dimensions.dims2indices), fu::typeof(DimensionalData.Dimensions.unalligned_dims2indices), lookups::Tuple{…}, dims::Tuple{…}, I::Tuple{…})
    @ DimensionalData.Dimensions ~/Documents/VectorDataCubes_Tutorial_EGU_2025/dev/DimensionalData/src/Dimensions/indexing.jl:91
 [20] split_alignments(fa::Function, fu::Function, dims::Tuple{…}, I::Tuple{…})
    @ DimensionalData.Dimensions ~/Documents/VectorDataCubes_Tutorial_EGU_2025/dev/DimensionalData/src/Dimensions/indexing.jl:89
 [21] dims2indices
    @ ~/Documents/VectorDataCubes_Tutorial_EGU_2025/dev/DimensionalData/src/Dimensions/indexing.jl:56 [inlined]
 [22] dims2indices
    @ ~/Documents/VectorDataCubes_Tutorial_EGU_2025/dev/DimensionalData/src/Dimensions/indexing.jl:32 [inlined]
 [23] _dim_getindex
    @ ~/Documents/VectorDataCubes_Tutorial_EGU_2025/dev/DimensionalData/src/array/indexing.jl:106 [inlined]
 [24] #getindex#176
    @ ~/Documents/VectorDataCubes_Tutorial_EGU_2025/dev/DimensionalData/src/array/indexing.jl:75 [inlined]
 [25] top-level scope
    @ REPL[69]:1
Some type information was truncated. Use `show(err)` to see complete types.

asinghvi17 · 2025-04-24T16:52:20Z

Fixed that.

I was thinking about feature collections a bit more, and I think we can do this:

If the input is a feature collection, decompose to a vector of features
If the input is a table and has a geometry column, then collect the output of Tables.rows into a vector of GeoInterface features
Implement GI.getgeom(geomlookup, i) which we can use everywhere

Then this is automatically able to subset as well. And getindex etc just work.

rafaqz · 2025-04-24T17:03:30Z

Sounds perfect.

The one downside is it will materialise a column table into rows, so allocates a lot. But users will learn to load the column if they only want the geoms

asinghvi17 · 2025-04-24T17:04:14Z

Yeah I don't see a way around rowtableizing unfortunately...maybe some lazy wrapper struct but that can DEFINITELY come later :D

rafaqz · 2025-04-24T17:08:46Z

Yeah had the EXACT same thought

rafaqz · 2025-04-24T17:09:50Z

We can also expand these back into columns in DimTable

asinghvi17 · 2025-04-24T17:14:04Z

Hm that would be interesting to do! Like an implicit dimstack

asinghvi17 · 2025-04-24T17:15:07Z

actually yeah, I like the dimstack idea more now that I think about it, it seems way more natural than this but also harder to do

felixcremer · 2025-04-25T12:51:56Z

Zonal on list of Tuples fails

Running zonal on the output of centroid I would have expected to get a Vector data cube with point entries.
Or is there a better way to get a point bsed vector data cube?

julia> zonal(identity, lai_da, of=Geometry(centroids))
ERROR: MethodError: no method matching npoint(::GeoInterface.PointTrait, ::Tuple{Float64, Float64})
The function `npoint` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  npoint(::GeoInterface.PentagonTrait, ::Any)
   @ GeoInterface ~/.julia/packages/GeoInterface/4tyIo/src/fallbacks.jl:83
  npoint(::GeoInterface.QuadTrait, ::Any)
   @ GeoInterface ~/.julia/packages/GeoInterface/4tyIo/src/fallbacks.jl:81
  npoint(::GeoInterface.TriangleTrait, ::Any)
   @ GeoInterface ~/.julia/packages/GeoInterface/4tyIo/src/fallbacks.jl:77
  ...

Stacktrace:
  [1] npoint(geom::Tuple{Float64, Float64})
    @ GeoInterface ~/.julia/packages/GeoInterface/4tyIo/src/interface.jl:175
  [2] _burn_geometry!(B::Raster{…}, ::GeoInterface.PointTrait, geom::Tuple{…}; shape::Nothing, verbose::Bool, boundary::Symbol, allocs::Rasters.Allocs{…}, fill::Bool, kw::@Kwargs{…})
    @ Rasters ~/Documents/VectorDataCubes_Tutorial_EGU_2025/dev/Rasters/src/methods/burning/geometry.jl:85
  [3] burn_geometry!(B::Raster{…}, data::Tuple{…}; kw::@Kwargs{…})
    @ Rasters ~/Documents/VectorDataCubes_Tutorial_EGU_2025/dev/Rasters/src/methods/burning/geometry.jl:4
  [4] boolmask!(dest::Raster{…}, data::Tuple{…}; invert::Bool, lock::Nothing, progress::Bool, threaded::Bool, burncheck_metadata::DimensionalData.Dimensions.Lookups.Metadata{…}, allocs::Rasters.Allocs{…}, geometrycolumn::Nothing, collapse::Rasters.NoKW, kw::@Kwargs{})
    @ Rasters ~/Documents/VectorDataCubes_Tutorial_EGU_2025/dev/Rasters/src/methods/mask.jl:354
  [5] boolmask(x::Tuple{Float64, Float64}; to::Tuple{X{…}, Y{…}}, invert::Bool, geometrycolumn::Nothing, kw::@Kwargs{})
    @ Rasters ~/Documents/VectorDataCubes_Tutorial_EGU_2025/dev/Rasters/src/methods/mask.jl:317
  [6] _mask(x::Raster{…}, with::Tuple{…}; kw::@Kwargs{})
    @ Rasters ~/Documents/VectorDataCubes_Tutorial_EGU_2025/dev/Rasters/src/methods/mask.jl:85
  [7] _mask(x::Raster{…}, with::Tuple{…})
    @ Rasters ~/Documents/VectorDataCubes_Tutorial_EGU_2025/dev/Rasters/src/methods/mask.jl:83
  [8] mask(x::Raster{…}; with::Tuple{…}, kw::@Kwargs{})
    @ Rasters ~/Documents/VectorDataCubes_Tutorial_EGU_2025/dev/Rasters/src/methods/mask.jl:76
  [9] _zonal(f::Function, x::Raster{…}, ::GeoInterface.PointTrait, geom::Tuple{…}; skipmissing::Bool, spatialslices::DimensionalData.Dimensions.Lookups._True, missingval::Missing, kw::@Kwargs{})
    @ Rasters ~/Documents/VectorDataCubes_Tutorial_EGU_2025/dev/Rasters/src/methods/zonal.jl:112
 [10] _zonal
    @ ~/Documents/VectorDataCubes_Tutorial_EGU_2025/dev/Rasters/src/methods/zonal.jl:103 [inlined]
 [11] _zonal
    @ ~/Documents/VectorDataCubes_Tutorial_EGU_2025/dev/Rasters/src/methods/zonal.jl:97 [inlined]
 [12] _alloc_zonal(f::Function, x::Raster{…}, geoms::DimVector{…}, n::Int64; spatialslices::DimensionalData.Dimensions.Lookups._True, missingval::Missing, kw::@Kwargs{…})
    @ Rasters ~/Documents/VectorDataCubes_Tutorial_EGU_2025/dev/Rasters/src/methods/zonal.jl:183
 [13] _zonal(f::Function, x::Raster{…}, ::Nothing, data::Geometry{…}; progress::Bool, threaded::Bool, geometrycolumn::Nothing, missingval::Missing, kw::@Kwargs{…})
    @ Rasters ~/Documents/VectorDataCubes_Tutorial_EGU_2025/dev/Rasters/src/methods/zonal.jl:141
 [14] _zonal
    @ ~/Documents/VectorDataCubes_Tutorial_EGU_2025/dev/Rasters/src/methods/zonal.jl:136 [inlined]
 [15] _zonal(f::Function, x::Raster{…}, of::Geometry{…}; kw::@Kwargs{…})
    @ Rasters ~/Documents/VectorDataCubes_Tutorial_EGU_2025/dev/Rasters/src/methods/zonal.jl:96
 [16] (::Rasters.var"#1047#1048"{…})(xo::Raster{…})
    @ Rasters ~/Documents/VectorDataCubes_Tutorial_EGU_2025/dev/Rasters/src/methods/zonal.jl:77
 [17] open(f::Rasters.var"#1047#1048"{…}, A::Raster{…}; kw::@Kwargs{})
    @ Rasters ~/Documents/VectorDataCubes_Tutorial_EGU_2025/dev/Rasters/src/array.jl:154
 [18] open(f::Function, A::Raster{…})
    @ Rasters ~/Documents/VectorDataCubes_Tutorial_EGU_2025/dev/Rasters/src/array.jl:149
 [19] #zonal#1046
    @ ~/Documents/VectorDataCubes_Tutorial_EGU_2025/dev/Rasters/src/methods/zonal.jl:76 [inlined]
 [20] top-level scope
    @ REPL[157]:1
Some type information was truncated. Use `show(err)` to see complete types.

asinghvi17 · 2025-04-25T13:09:20Z

Extract will do it but it won't give you a VDC. Guess we have to get that next :)

But zonal probably doesn't support point geoms because it tries to burn them.

rafaqz · 2025-04-25T13:19:47Z

Yeah extract is better for points, zonal doesn't have the point optimisations

asinghvi17 · 2025-04-26T19:39:56Z

@rafaqz extract is now also multidimensional (but VERY inefficient)!

asinghvi17 added 5 commits January 8, 2025 10:26

WIP: geometry lookup

5dcb1f7

Multi/polygon writing

4bccf40

Add crs, reproject

6717902

Hook up CRS, export GeometryLookup

8d6f670

add super basic tests

b614376

Merge branch 'main' into as/vectordatacubes

5ec0af6

asinghvi17 added 3 commits April 22, 2025 11:23

move geometry_lookup to a folder

ab7d427

add overrides for zonal when it's passed a GeometryLookup or dim that…

1c55155

… wraps that

make zonal use the fast-path missing return only if skipmissing = true

55cb66d

This way we can force skipmissing=false and manually do it, if we want e.g. vector data cubes.

asinghvi17 commented Apr 22, 2025

View reviewed changes

implement automatic spatial slicing for datacubes

5f6b368

with style and ease

asinghvi17 added 2 commits April 22, 2025 15:39

add comments, set spatialslices to _True initially for assured const

09d8ac3

fix geomlookup method

3757a35

asinghvi17 added 11 commits April 22, 2025 15:51

wip fancy rebuild / dim detection

aa1bcc5

implement isfalse

debb290

so that we can check that the thing is explicitly false, this is quite useful for spatialslice because you can provide anything but false there, and it assumes true

add a missingval kwarg so we can fix things correctly

16ef85b

missingval support

af2e447

fix run: remove static, let threaded be anything for extensibility

233018f

pass spatial slices through the geometrylookup zonal

41ef26e

disallow missings in geom lookup

1b7d26c

allow metadata kwarg to rebuild

19775da

we need an interface / spec for what rebuild must accept in DD!!

plotting support for dim 1 vectordatacubes

97474fe

default missingval is missing not nothing

e1ab26d

Geometry named dim

b0293ca

asinghvi17 added 2 commits April 24, 2025 06:09

fix

3bd7451

make crs better on rasters

f89e23b

asinghvi17 added 2 commits April 24, 2025 10:12

allow Near on points

2e44421

enable at on points / forward for other things

2184b92

asinghvi17 added 4 commits April 26, 2025 10:31

fix touches

88230bb

minor fixes

626b8ee

make extract multidimensional

899b36a

make the change non breaking

26ce1a0

asinghvi17 added 10 commits April 26, 2025 16:07

Clean up lookups

a97ae9f

Update methods.jl

056cd0d

Update zonal.jl

3b052b9

only run zonal if not on fast, all-is-missing path

629cf6a

basic netcdf io, could be more optimized / streamlike

4a56b50

fix

3c1cc29

add some test files as text that can be fed into ncgen

71700f5

add decoders for all geomtraits

eb6fa75

use the test files in nc_io.jl which should now be a test file

e9f5186

fully implement, and "test" all six geom types

e94b28b

WIP: geometry lookup #854

Are you sure you want to change the base?

WIP: geometry lookup #854

Uh oh!

Conversation

asinghvi17 commented Jan 8, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

felixcremer commented Apr 22, 2025

Uh oh!

asinghvi17 commented Apr 22, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

rafaqz commented Apr 22, 2025

Uh oh!

asinghvi17 Apr 22, 2025

Choose a reason for hiding this comment

Uh oh!

asinghvi17 Apr 22, 2025

Choose a reason for hiding this comment

Uh oh!

asinghvi17 commented Apr 22, 2025

Uh oh!

rafaqz commented Apr 22, 2025

Uh oh!

asinghvi17 commented Apr 22, 2025

Uh oh!

felixcremer commented Apr 24, 2025

Uh oh!

asinghvi17 commented Apr 24, 2025

Uh oh!

rafaqz commented Apr 24, 2025

Uh oh!

asinghvi17 commented Apr 24, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

rafaqz commented Apr 24, 2025

Uh oh!

rafaqz commented Apr 24, 2025

Uh oh!

asinghvi17 commented Apr 24, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

asinghvi17 commented Apr 24, 2025

Uh oh!

felixcremer commented Apr 25, 2025

Zonal on list of Tuples fails

Uh oh!

asinghvi17 commented Apr 25, 2025

Uh oh!

rafaqz commented Apr 25, 2025

Uh oh!

asinghvi17 commented Apr 26, 2025

Uh oh!

Uh oh!

asinghvi17 commented Jan 8, 2025 •

edited

Loading

asinghvi17 commented Apr 22, 2025 •

edited

Loading

asinghvi17 commented Apr 24, 2025 •

edited

Loading

asinghvi17 commented Apr 24, 2025 •

edited

Loading