single point extraction with COPERNICUS_30 returns incorrect results

[JeroenVerstraelen]

The following simple example shows that single point extraction for COPERNICUS_30 always produces "[[0,0]]" as output:
https://gist.github.com/JeroenVerstraelen/327362be431b8158fcabb7b5d10f3d66

The same code does work for the SENTINEL2_L2A collection.

[JeroenVerstraelen]

Note that it gives the correct output when using 2 point features in a FeatureCollection, but it returns 0.0 if you use 1 point feature instead.

[JeroenVerstraelen]

# Multiple points:
  # combinedRDD == (SpatialKey(0,0),(ArrayMultibandTile(128,128,1,float32),CompactBuffer(Feature(POINT (10.8774075137 46.7888937363),0), Feature(POINT (10.8674075137 46.7988937363),1))))
  # datacube.metadata = TileLayerMetadata(float32,LayoutDefinition(Extent(10.867317770256486, 46.7278723686927, 10.938428881367598, 46.79898347980381),CellSize(2.7777777777777957E-4,2.777777777777657E-4),2x2 tiles,256x256 pixels),Extent(10.8673177703, 46.7888039928, 10.8774972572, 46.7989834798),EPSG:4326,KeyBounds(SpatialKey(0,0),SpatialKey(0,0)))
# Single point:
  # combinedRDD == (SpatialKey(0,0),(ArrayMultibandTile(32,32,1,float32),CompactBuffer(Feature(POINT (10.8774075137 46.7888937363),0))))
  # datacube.metadata == TileLayerMetadata(float32,LayoutDefinition(Extent(10.877317770256486, 46.78019288168888, 10.886206659145374, 46.78908177057777),CellSize(2.777777777777657E-4,2.777777777778212E-4),1x1 tiles,32x32 pixels),Extent(10.8773177703, 46.7888039928, 10.8774972572, 46.7889834798),EPSG:4326,KeyBounds(SpatialKey(0,0),SpatialKey(0,0)))

Single point:
datacube.collect().head._2.band(0) = PaddedTile(ArrayTile(1,1,float32),0,0,32,32)
datacube.collect().head._2.band(0).toArrayDouble() == 
0 = 0.0
1 = NaN
2 = NaN
...
1024 = Nan

The singular value in the ArrayTile is 0.0, the rest is padded with NaN values.

So the issue appears to occur when setting up the extent for the entire datacube.

[JeroenVerstraelen]

_load_collection_cached::spatial_extent == {'crs': 'EPSG:4326', 'east': 10.877497257232971, 'north': 46.78898347980381, 'south': 46.78880399282733, 'west': 10.877317770256486}
PyramidFactory.datacube_seq::polygons == MULTIPOLYGON (((10.8774972572 46.7888937363, 10.877470972 46.7888302781, 10.8774075137 46.7888039928, 10.8773440555 46.7888302781, 10.8773177703 46.7888937363, 10.8773440555 46.7889571945, 10.8774075137 46.7889834798, 10.877470972 46.7889571945, 10.8774972572 46.7888937363)))
PyramidFactory.datacube_seq::polygons.areaInSquareMeters == 193.35853276962737

readKeysToRasterSourcesResult :: GeoTiffReprojectRasterSource(/eodata/auxdata/CopDEM/COP-DEM_GLO-30-DGED/DEM1_SAR_DGE_30_20110517T170701_20140817T170857_ADS_000000_bma2.DEM/Copernicus_DSM_10_N46_00_E010_00/DEM/Copernicus_DSM_10_N46_00_E010_00_DEM.tif, EPSG:4326, TargetRegion(GridExtent(Extent(10.877317770256486, 46.7888039928, 10.877595548034263, 46.78908177057777), CellSize(2.777777777777657E-4,2.777777777778212E-4), 1x1)), NearestNeighbor)

The rastersource has the following gridExtent:
gridExtent = {geotrellis.raster.GridExtent@14626} GridExtent(Extent(10.877317770256486, 46.7888039928, 10.877595548034263, 46.78908177057777), CellSize(2.777777777777657E-4,2.777777777778212E-4), 1x1)
 extent: geotrellis.vector.Extent  = {geotrellis.vector.Extent@14637} Extent(10.877317770256486, 46.7888039928, 10.877595548034263, 46.78908177057777)
 cellwidth: double  = 2.777777777777657E-4
 cellheight: double  = 2.777777777778212E-4
 cols: java.lang.Object  = {java.lang.Long@14638} 1
 rows: java.lang.Object  = {java.lang.Long@14638} 1

This happens because in readKeysToRasterSources:
reprojectedBoundingBox.extent == Extent(10.8773177703, 46.7888039928, 10.8774972572, 46.7889834798)
reprojectedBoundingBox.extent.area == 3.221556522110769E-8
alignedExtent == GridExtent(Extent(10.877317770256486, 46.7888039928, 10.877595548034263, 46.78908177057777), CellSize(2.777777777777657E-4,2.777777777778212E-4), 1x1) # 1 col, 1 row
alignedExtent.extent == Extent(10.877317770256486, 46.7888039928, 10.877595548034263, 46.78908177057777)
alignedExtent.extent.area == 7.716049382619775E-8

dry_run

0 = {DataTrace} <DataTrace#125825649973264(#125825649934064, aggregate_spatial, {'geometries': DriverVectorCube(\ngeometries:\n                 ...tor_cube_dummy:  True\ncrs: 'EPSG:4326'\nbbox: (10.877407513744728, 46.78889373631557, 10.877407513744728, 46.78889373631557))})>
 children = {list: 0} []
 parent = {DataTrace} <DataTrace#125825649934064(#125825687748496, weak_spatial_extent, {'west': 10.877317770256486, 'south': 46.78880399282733, 'east': 10.877497257232971, 'north': 46.78898347980381, 'crs': 'EPSG:4326'})>

The issue arises because of the hardcoded 10m point buffering when calculating the weak_spatial_extent in the dry_run.
An issue for this exists already here:
Open-EO/openeo-python-driver#148

[JeroenVerstraelen]

Summary:

• This issue only occurs when using a single point as weak_spatial_extent for layers with a resolution lower than 10m
• A proper fix is to replace the hardcoded 10m point buffering with the actual resolution of the datacube
• A workaround is to provide more than one point (as long as both points do not fall within the same pixel)
• Another workaround is to provide a spatial_extent in load_collection after this PR is merged: fix extent order when calculating global_extent #358 openeo-python-driver#359

[jdries]

Note that we want to avoid telling users to use a spatial extent in load_collection. That would be opposite to other advice we're giving, causing confusion.