Investigation and clarification needed on /cog/statistics endpoint for area of interest (AOI). Data statistics do not appear consistent.

[kyle-lesinger]

Issue: Maximum value is outside of range using /cog/statistics for an area of interest (AOI). (subset of the COG file)

Here is some information in potential errors or uncertainties when using /cog/statistics endpoint within raster API.

This is for the notebook user_data_notebooks/cms-global-map-mangrove_User_Notebook.ipynb. For my specific notebook with mangrove data for GHG, I've noticed that the maximum value within the selected AOI is outside of the bounds of the maximum value in the tile. This makes me think that /cog/statistics is selecting some other data, or I may just be mis-interpreting the API functionality so I would need some clarification.

Here is the PR in ghgc-docs : mangrove-data.
This is the STAC collection

1.) For AOI selected within the Australia tile, the maximum value in the collection statistics is 145.98 (as seen in the collection id cms-mangrove-agb-canopyheight-grid-v1.3-Australia for variable agb). But running the generate_stats function which calls POST /cog/statistics endpoint, the maximum value is 212.61 which is outside of the range. Even though its a subset of the tile, it should not be >145.98.

2.) I looked into all 3 assets (agb, hmax95, hba) to see if it was getting those statistics, but none of them have values >145.98.

3a.) I looked briefly through other User_notebooks and I believe that I found another inconsistency within odiac-ffco2-monthgrid-v2022_User_Notebook.ipynb in ghgc-docs. Here is the STAC collection/items link. Within this notebook for item_id odiac-ffco2-monthgrid-v2022-202112, the maximum value from generate_stats is 404594. But the maximum value within the item['statistics'] in STAC collection is 2497 for the asset co2-emissions.

3b.) To slightly complicate things, when I looked at the .tif file for 202112 in COG Viewer, I see that the maximum is around 585 (not at 404594 or 2497). Here is the file path s3://ghgc-data-store/odiac-ffco2-monthgrid-v2022/odiac2022_1km_excl_intl_202112.tif. The histogram bins may influence what we perceive as the maximum but after changing the rescale there are no values >100000.

This may require more information, so please reach out to me and I'll offer any additional assistance or insight you may need.

[vincentsarago]

@kyle-lesinger Can you share how you're getting the statistics from the /cog/statistics endpoint

In your notebook I can see that you're getting the statistics from the high resolution data while in the /cog/statistics endpoint the API will mostly use the overview within the COGs. Depending how the overviews were created (resampling) the overviews might have slightly different values.

you can easily get the stats using rio-tiler such as:

from rio_tiler.io import Reader

with Reader("cog.tif") as src:
    stats = src.statistics()
    print(stats["b1"].max)

    # Force reading the stats from the high resolution
    stats_native = src.statistics(max_size=None)
    print(stats_native["b1"].max)

To slightly complicate things, when I looked at the .tif file for 202112 in COG Viewer, I see that the maximum is around 585 (not at 404594 or 2497). Here is the file path s3://ghgc-data-store/odiac-ffco2-monthgrid-v2022/odiac2022_1km_excl_intl_202112.tif. The histogram bins may influence what we perceive as the maximum but after changing the rescale there are no values >100000.

👇 explain the difference of statistics based on overview level

# File info
curl https://earth.gov/ghgcenter/api/raster/cog/info\?url\=s3://ghgc-data-store/odiac-ffco2-monthgrid-v2022/odiac2022_1km_excl_intl_202112.tif | jq      
{
  "bounds": [
    -180.0,
    -90.0,
    180.0,
    90.0
  ],
  "minzoom": 0,
  "maxzoom": 4,
  "band_metadata": [
    [
      "b1",
      {}
    ]
  ],
  "band_descriptions": [
    [
      "b1",
      "band_data"
    ]
  ],
  "dtype": "float32",
  "nodata_type": "Nodata",
  "colorinterp": [
    "gray"
  ],
  "driver": "GTiff",
  "count": 1,
  "width": 43200,
  "height": 21600,
  "overviews": [
    2,
    4,
    8,
    16,
    32,
    64,
    128
  ],
  "nodata_value": -9999.0
}

# Max Size for each Overview
overview_factors = [
    2,
    4,
    8,
    16,
    32,
    64,
    128
]
# max size `max(height, width)`
size = 43200
ii = 0
for level in overview_factors:
    print(f"Overview Level {len(overview_factors) - ii - 1} (dec={level})", size // leve)
    ii += 1

>>> Overview Level 6 (dec=2) 21600
Overview Level 5 (dec=4) 10800
Overview Level 4 (dec=8) 5400
Overview Level 3 (dec=16) 2700
Overview Level 2 (dec=32) 1350
Overview Level 1 (dec=64) 675
Overview Level 0 (dec=128) 337

# Default Statistics (max_size=1024) - Will read overview level 2 (decimation=32)
curl https://earth.gov/ghgcenter/api/raster/cog/statistics\?url\=s3://ghgc-data-store/odiac-ffco2-monthgrid-v2022/odiac2022_1km_excl_intl_202112.tif | jq '.b1.max'
2497.01904296875

# Statistics max_size=2048 - Will read overview level 3 (decimation=16)
curl https://earth.gov/ghgcenter/api/raster/cog/statistics\?url\=s3://ghgc-data-store/odiac-ffco2-monthgrid-v2022/odiac2022_1km_excl_intl_202112.tif\&max_size\=2048 | jq '.b1.max'
5205.00732421875

# Statistics max_size=512 - Will read overview level 1 (decimation=64)
curl https://earth.gov/ghgcenter/api/raster/cog/statistics\?url\=s3://ghgc-data-store/odiac-ffco2-monthgrid-v2022/odiac2022_1km_excl_intl_202112.tif\&max_size\=512 | jq '.b1.max'
907.970458984375

Preview (max_size=1024)

url: https://earth.gov/ghgcenter/api/raster/cog/preview?url=s3://ghgc-data-store/odiac-ffco2-monthgrid-v2022/odiac2022_1km_excl_intl_202112.tif&rescale=0,100&max_size=1024&colormap_name=viridis

Preview (max_size=512)

url: https://earth.gov/ghgcenter/api/raster/cog/preview?url=s3://ghgc-data-store/odiac-ffco2-monthgrid-v2022/odiac2022_1km_excl_intl_202112.tif&rescale=0,100&max_size=512&colormap_name=viridis

Looking at the data at different overview level, I'm pretty sure the overviews were created using non-nearest resampling but maybe a sum algorithm, which will explain the statistics difference at each level

[kyle-lesinger]

@kyle-lesinger Can you share how you're getting the statistics from the /cog/statistics endpoint

In your notebook I can see that you're getting the statistics from the high resolution data while in the /cog/statistics endpoint the API will mostly use the overview within the COGs. Depending how the overviews were created (resampling) the overviews might have slightly different values.

Sorry for the delay. I'm getting the statistics for the mangrove dataset within ghgc-docs/ghgc_utils.py. This script will

• post a request using the area of interest (as the geojson), the s3 url where the .tif is located

    # A POST request is made to submit the data associated with the item of interest (specific observation) within the boundaries of the polygon to compute its statistics
    result = requests.post(

        # Raster API Endpoint for computing statistics
        f"{kwargs['url']}/cog/statistics",

        # Pass the URL to the item, asset name, and raster identifier as parameters
        params=params,   -- this is the s3 url of the .tif file

        # Send the GeoJSON object (polygon) along with the request
        json=geojson, ---- area of interest which contains the bounding box specified within the notebook

    # Return the response in JSON format
    ).json() 

Then the function further processes the data to return a pandas dataframe with the statistics

datetime | min | max | mean | count | 
2000-01-01T00:00:00Z | 57.41016006469726562500 | 212.61692810058593750000 | 121.31508636474609375000 | 618041.00000000000000000000

[kyle-lesinger]

👇 explain the difference of statistics based on overview level

# File info
curl https://earth.gov/ghgcenter/api/raster/cog/info\?url\=s3://ghgc-data-store/odiac-ffco2-monthgrid-v2022/odiac2022_1km_excl_intl_202112.tif | jq      
{
  "bounds": [
    -180.0,
    -90.0,
    180.0,
    90.0
  ],
  "minzoom": 0,
  "maxzoom": 4,
  "band_metadata": [
    [
      "b1",
      {}
    ]
  ],
  "band_descriptions": [
    [
      "b1",
      "band_data"
    ]
  ],
  "dtype": "float32",
  "nodata_type": "Nodata",
  "colorinterp": [
    "gray"
  ],
  "driver": "GTiff",
  "count": 1,
  "width": 43200,
  "height": 21600,
  "overviews": [
    2,
    4,
    8,
    16,
    32,
    64,
    128
  ],
  "nodata_value": -9999.0
}

# Max Size for each Overview
overview_factors = [
    2,
    4,
    8,
    16,
    32,
    64,
    128
]
# max size `max(height, width)`
size = 43200
ii = 0
for level in overview_factors:
    print(f"Overview Level {len(overview_factors) - ii - 1} (dec={level})", size // leve)
    ii += 1

>>> Overview Level 6 (dec=2) 21600
Overview Level 5 (dec=4) 10800
Overview Level 4 (dec=8) 5400
Overview Level 3 (dec=16) 2700
Overview Level 2 (dec=32) 1350
Overview Level 1 (dec=64) 675
Overview Level 0 (dec=128) 337

# Default Statistics (max_size=1024) - Will read overview level 2 (decimation=32)
curl https://earth.gov/ghgcenter/api/raster/cog/statistics\?url\=s3://ghgc-data-store/odiac-ffco2-monthgrid-v2022/odiac2022_1km_excl_intl_202112.tif | jq '.b1.max'
2497.01904296875

# Statistics max_size=2048 - Will read overview level 3 (decimation=16)
curl https://earth.gov/ghgcenter/api/raster/cog/statistics\?url\=s3://ghgc-data-store/odiac-ffco2-monthgrid-v2022/odiac2022_1km_excl_intl_202112.tif\&max_size\=2048 | jq '.b1.max'
5205.00732421875

# Statistics max_size=512 - Will read overview level 1 (decimation=64)
curl https://earth.gov/ghgcenter/api/raster/cog/statistics\?url\=s3://ghgc-data-store/odiac-ffco2-monthgrid-v2022/odiac2022_1km_excl_intl_202112.tif\&max_size\=512 | jq '.b1.max'
907.970458984375

Thank you for providing a detailed description of how differing zoom levels which cause difference in the statistics reported, this makes sense.

Let me think about it over the weekend because when individuals are physically pulling data within the notebook or they choose to do analysis with the data from our STAC/buckets, then this discrepancy between what is posted on the legend and what is reported in legend will cause some uncertainty.

Thanks again for the info.

[kyle-lesinger]

@siddharth0248 @slesaad @Jeanne-le-Roux @paridhi-parajuli See message above regarding differences in /cog/statistics values based on decimation. There is nothing wrong with /cog/statistics but there is a discrepancy sometimes between what the legend shows and what the actual value is.

I can reproduce similar stats that were being obatined by the generates_stats function by calling different decimation sizes. We can discuss this at another time. Maybe one solution would be in notebooks to call a specific level to generate the highest resolution stats when plotting the color bar?

> curl https://earth.gov/ghgcenter/api/raster/cog/statistics\?url\=s3://ghgc-data-store/cms-mangrove-agb-canopyheight-grid-v1.3/Mangrove_agb_Australia_2000year.tif\&max_size\=1024 | jq '.b1.max'
145.98110961914062

> curl https://earth.gov/ghgcenter/api/raster/cog/statistics\?url\=s3://ghgc-data-store/cms-mangrove-agb-canopyheight-grid-v1.3/Mangrove_agb_Australia_2000year.tif\&max_size\=2048 | jq '.b1.max'
166.1106414794922

> curl https://earth.gov/ghgcenter/api/raster/cog/statistics\?url\=s3://ghgc-data-store/cms-mangrove-agb-canopyheight-grid-v1.3/Mangrove_agb_Australia_2000year.tif\&max_size\=4096 | jq '.b1.max'
214.48577880859375

here is the image in question

And here is the maximum value reported for that area of interest 212.61692810058593750000