Simple KD tree is insufficient for coordinate transformation

[danlooo]

We want fast transformation between cell ids and lan/lot coordinates.
Since the DGGRID polygons are regular, we can just store the center points in a KD-Tree and assume a voronoi like partition.
However, this seems to be only true at the planar face of the DGGS polyhedron as demonstrated in this R script using centers and boundaries exported from DGGRID as GeoJSON:

library(tidyverse)
library(sf)

preset <- "isea4h"
centers <- str_glue("centers/{preset}.geojson") |> st_read() |> st_as_sf()
boundaries <- str_glue("boundaries/{preset}.geojson") |> st_read() |> st_as_sf()
voronoi <- centers |> st_combine() |> st_voronoi()

ggplot() +
  geom_sf(data = centers) +
  geom_sf(data = boundaries, aes(color = "Boundaries"), fill = NA,) +
  geom_sf(data = voronoi, aes(color = "Voronoi"), fill = NA) +
  coord_sf(xlim = c(-60,60), ylim = c(0,90)) +
  scale_color_manual(values = c("Boundaries" = "red", "Voronoi" = "blue"))

The difference, of course, is the highest near the poles using high shape distorted ISEA projection.

In Addition, the cell_id of some points are different if they are calculated based on DGGRID boundary polygons or voronoi like partitioning of getting the nearest center point neighbor using R function nngeo::st_nn with metric sf::st_distance:

expand_grid(
  x = seq(-180,180, 10),
  y = seq(-90,90, 10)
) |>
  st_as_sf(crs = 4326, coords = c("x", "y")) |>
  st_join(boundaries |> rename(boundary_cell_id = cell_id), join = st_within) |>
  st_join(centers |> rename(centers_cell_id = cell_id), join = nngeo::st_nn) |>
  count(boundary_cell_id == centers_cell_id)
# A tibble: 2 × 3
#  `boundary_cell_id == centers_cell_id`     n                                                                                  geometry
#* <lgl>                                 <int>                                                                          <MULTIPOINT [°]>
#1 FALSE                                    20 ((60 -20), (90 0), (60 0), (60 20), (-20 50), (-80 50), (-80 80), (160 50), (100 80), ...
#2 TRUE                                    683 ((-40 -30), (-40 -20), (-30 -30), (-30 -20), (-30 -40), (-20 -40), (-10 -40), (-10 -30...

Do we really need system calls to DGGRID binary for a single point to cell id look up?
Seems we need to re-project not to the entire sphere, but only to the corresponding face of the DGGS polyhedron.
This would involve to create a custom metric type of Distances.jl that is doing the re-projection to the corresponding polyhedron face.
Am I right, @allixender ?

[danlooo]

DGGRID offers conversion between geographical coordinates and x,y on the polyhedron faces (PROJTRI). Two grid points for each of the 20 faces of the DGGS icosahedron and their corresponding voronoi partition:

Points are equally spaced in lon/lat space:

Cell center points:

Seems that we can still build a KD-Tree in the projtri space.

[allixender]

Hi @danlooo , just a comment, you are doing interesting stuff. Generally I would want that we can do single API calls into DGGRID for single point conversions, just like H3 and S2 are offering it. The thing of course with DGGRID as well as with e.g. rHEALPix or OpenEAGGR is that the correct DGGS model needs to be instantiated first before you can ask it to do any calculations. H3 and S2 have a fixed rotation and and single geometry so single point conversions are always operating already on the hard-coded DGGS model.

[danlooo]

Just to make sure, boundary points projected by DGGRID and voronoi partition of center points are identical in PROJITRI space, because this is before the projection where the cells are regular.

[asinghvi17]

Purely out of curiosity, did you look into working directly in Cartesian coordinates (x, y, z) after transforming lat/long to e.g a unit sphere for this? It seems to me that you could then avoid any projections at all and work in 3D space on the surface of the sphere, where all polygons would have to be more or less regular.

[danlooo]

@asinghvi17 Indeed, there is no difference in terms of accuracy in using either polar or Cartesian coordinates. It's all about finding an efficient data structure for billions or even trillions of cells (We target high-res satellite raster data here).
kD trees are inaccurate (works only in the plane) and inefficient (tree and not tensor, we treat raster data as vector data here).
Cartesian data cubes of a sphere surface are very sparse and thus also inefficient.
This is why DGGS.jl currently uses DGGRID Q2DI index instead.