how xr_quantile computes percentiles with overlap from different times and NaNs

[ithaca-oyler]

Hi again! This is the continued discussion from opendatacube/odc-stac#177. As mentioned in the other issue, I'm using HLS data (both HLSL30.v2.0 and HLSS30.v2.0 collections) for a short time frame (5-16-2024 through 5-18-2024) over a small region of interest in Washington state. I am currently using odc-geo==0.4.6, odc-stac==0.3.9, and odc-algo=0.2.2. We're generating composites for the 15th, 50th, and 90th percentiles. I've attached a jupyter notebook with a reproducible example.

As mentioned in the odc-stac issue, using a slightly different method, we anticipate the percentile composites to look something like this (the 50th percentile is shown here).

However, in using xr_quantile, we're getting different outputs. This is what the RGB of cloud masked input data looks like, one time slice per day.

This is what the 15th, 50th, and 90th percentiles look like when generated with xr_quantile. Given how much coverage there is with the input imagery over the three days, we would expect each percentile image to have lots of coverage with at least one valid value per pixel (as opposed to the widespread NaN values, which are mostly present in the 90th percentile).

Based on these results, I have several questions about the algorithm. Starting with a simple example, if there are 3 scenes that overlap for a particular pixel, each from a different day, a pixel could have up to 3 values. Let's say those three values are 100, 500, and 1000. With the way we compute percentiles, the 15th, 50th, and 90th percentiles would be 220, 500, and 900, respectively.

1. How does xr_quantile handle NaN values in the pixel value arrays? For the three day example, (nan, 500, 1000), np.percentile() would return all NaN, but np.nanpercentile() would return 15th=575, 50th=750, 90th=950. How would xr_quantile handle the (nan, 500, 1000) array?

2. Related to the first question, is this value what the nodata argument in xr_quantile is designed as an input for? For example, if I use a nodata=0, the array (nan, 500, 1000) would become (0, 500, 1000), and the outputs would be 15th=150, 50th=500, 90th=900.

3. Do we have any idea what might be going on with the overlap at the different quantiles? As I mentioned, we would expect all three days worth of pixel values to be accounted for in each quantile calculation so that each image has coverage like quantile=0.15 does.

replication-code-odc-algo.zip
Question -- what's the best way to provide you with the sample dataset? It's too big to attach here.

[maawoo]

Have you tried calculating the quantiles using Xarray's methods instead of the implementation of this package? Just make sure to also install the package Numbagg in order to speed things up. Otherwise you will run into issues that have been around for some while and probably led to this package's implementation in the first place...

[robbibt]

Hey @ithaca-oyler, I'll try and have a look at this this week! Just to clarify, have you seen different results with respect to NaN handling from xr_quantile vs. the native .quantile() functionality from xarray? It would be great to do a side-by-side comparison to verify that the results you're seeing are different from that existing functionality (although .quantile() is extremely slow for large datasets given it loops through each pixel one-by-one).

(Just pinging the original algorithm developer @kieranricardo as well as @omad and @alexgleith as I believe they have been working with this code base recently too and may know if there's a more recent implementation that might avoid some of these issues)

[alexgleith]

Hey folks

There were a few nan bugs in the geomedian code, which I think have been resolved. That's where I'd look. Maybe it's sorting including nan, and nan is up the top, so the 90th percentile is nan.

I guess dropping nans and taking a percentile of that is the right approach, unless they're all nan... not sure though.

Some ref to changes: auspatious/datacube-compute@4e03258#diff-f5e0ecb72dacb41ee8b0a629b23efb91d1a88c95b2c39272b48ade4830ec0938L61

[ithaca-oyler]

Hi everyone! Thanks for taking a look at this issue. I have some more insights about my issue in particular.

It turns out by looking at three days, I made an edge case! I plotted the valid_counts data array created in np_percentile() and as shown here, the yellow stripe shows where there is an overlap in data across all three days. With the valid_counts input into keep_good_np() set to >=3, we were only getting outputs in the yellow stripe.

By setting valid_counts >=0, I got the expected outputs (the 15th percentile is shown here). I confirmed there is indeed a difference between the rasters by doing a raster subtraction between the 90th and 15th percentiles (the second image shown here). Where the image is brighter is also where the input imagery overlaps the most. For my use case, this likely won't be a big issue at scale, but it is an issue for small time frames and in areas with sparse imagery (<=3 scenes).

Now, curiously enough, the issue persists when the nodata value in xr_quantile is set to 0. I can understand how it would result in lower quantile values, but not 0s where valid_counts < 3.

In response to the comments above:

• @maawoo Thanks for pointing me to Numbagg and the Xarray methods! I'll have to give them a shot to see if it speeds things up or changes things
• @robbibt I haven't checked the NaN handling between xr_quantile and .quantile() with xarray yet, but I'll be sure to update you when I do the comparison.
• @alexgleith Do you think this is also what could be happening when the nodata value is set to 0 in xr_quantile?