WIP: geometry lookup

[asinghvi17]

TODOs:

• CRS support
  • Make reproject reproject the geometry lookups too
  • In-memory resample using e.g zonal, to rasterize a VDC
• Allow feature attributes to persist within geometry lookups
  • candidate 1: maybe have a metadata field?
  • candidate 2: table like feature collection structure, with metadata being a namedtuple of vectors...
• Plot recipes - plottype on a dimvector with geometry lookup should check the trait and return poly. Then poly should have a recipe on a dimvector with geometry lookup.
• Make everything return vectors, except when you ask for points
• Add I/O support via CF conventions
  • Read: point, line, polygon
    • Actual reading given a dict
    • Hook up to CF
  • Write: point, line, polygon
    • [] Actual writing given a dict
    • Hook up to CF
• Since we're depending on GeometryOps, use GeometryOps functions internally
• Add test cases based on xvec tests
• Translate all xvec examples

[felixcremer]

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

[asinghvi17]

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]

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

[asinghvi17]

@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.

[asinghvi17]

from Slack this is extremely cool and we should do this

[asinghvi17]

I guess the place to do IO is:

• _cdmlookup for read
• _def_dim_var! for write
  ?

[rafaqz]

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

[asinghvi17]

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

[rafaqz]

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]

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

[rafaqz]

Yeah had the EXACT same thought

[rafaqz]

We can also expand these back into columns in DimTable

[asinghvi17]

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

[asinghvi17]

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]

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]

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]

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

[asinghvi17]

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

[lazarusA]

any chance to merge things that work, and then leave complex ones for a different PR? Basic functionally for vectordatacubes at least would be nice to support.

[asinghvi17]

This might be breaking for zonal I think. But would be nice to have esp. leading into SDSL.

[asinghvi17]

iirc this also needs rafaqz/DimensionalData.jl#991

[lazarusA]

Wow, big changes everywhere 😆. I'm testing this branch again, let's see. Happy to test other changes 😀

[asinghvi17]

Ok that dd pr was at least merged, let me get this resolving and maybe add some ci

[asinghvi17]

But @rafaqz is unlikely to be able to help much for the next 2 months or so. Maybe we can just push this through if he's ok with it and make sure it is not breaking unless you have opted into geometrylookup somehow? lock it behind some flag or something?