Spatial partitioning, sorting and shuffling

When dealing with large sets of geometries it would be nice if we could partition (chunk) the geometry coordinate and the GeometryIndex based on spatial locality (thus requiring spatial sorting or shuffling), like explained for geo-dataframes in dask-geopandas' spatial partitioning user guide.

This would require a good amount of work both here and upstream, though:

• refactor xarray.IndexVariable to allow dask arrays or other lazy arrays https://github.com/pydata/xarray/pull/8124 (although we can work around this by using regular xarray.Variable objects for now)
• add some sort of generic, out-of-core index in Xarray https://github.com/pydata/xarray/issues/1650, or add our own ad-hoc version here
• compute space-filling curve values from geometries (copied or factored out from dask-geopandas?)
• implement spatial shuffling, probably reusing what is being added in Xarray https://github.com/pydata/xarray/pull/9320

Start with the dask-geopandas spatial partitioning guide, then review the linked xarray.IndexVariable PR and out-of-core index issue. Check the linked Xarray spatial-shuffling PR and determine how space-filling curve values and chunked geometry indexes should fit in xvec. Done means spatially partitioned geometry coordinates and GeometryIndex support the intended out-of-core sorting and shuffling.