Create a dataset for benchmarking¶

Combine three datasets:

1. The processed OCR burn probability data
2. The original burn probability data from Riley et al. (to create a "non-burnable" mask)
3. Historical fire perimeters

Note, this script has not been optimized for performance. It requires at least 400GB of memory (I used m8g.48xlarge)

In [1]:

import gc
import os

import geopandas as gpd
import icechunk
import numpy as np
import rasterio.features
import xarray as xr

# import seaborn as sns    # plotting
from ocr import catalog  # for the riley data

import gc import os import geopandas as gpd import icechunk import numpy as np import rasterio.features import xarray as xr # import seaborn as sns # plotting from ocr import catalog # for the riley data

Read in and pre-process the data¶

OCR burn probability¶

In [2]:

version = '0.8.0'
setup = 'production'

storage = icechunk.s3_storage(
    bucket='carbonplan-ocr',
    prefix=f'output/fire-risk/tensor/{setup}/v{version}/ocr.icechunk',
    from_env=True,
)
repo = icechunk.Repository.open(storage)
session = repo.readonly_session('main')

ds = xr.open_zarr(session.store, consolidated=False, zarr_format=3)

version = '0.8.0' setup = 'production' storage = icechunk.s3_storage( bucket='carbonplan-ocr', prefix=f'output/fire-risk/tensor/{setup}/v{version}/ocr.icechunk', from_env=True, ) repo = icechunk.Repository.open(storage) session = repo.readonly_session('main') ds = xr.open_zarr(session.store, consolidated=False, zarr_format=3)

Riley et al. burn probability¶

In [3]:

# --- read in
riley_2011_30m_4326 = catalog.get_dataset('2011-climate-run-30m-4326').to_xarray()[['BP']]
# there are slight mismatches with the coordinates in `ds`
# interpolate riley data to exactly match ds coordinates
riley_interp = riley_2011_30m_4326.interp(
    latitude=ds.latitude, longitude=ds.longitude, method='nearest'
)
# assign coordinates
ds['riley_BP_2011'] = riley_interp.BP
# create burnable mask
ds['riley_burnable_mask'] = xr.where(ds['riley_BP_2011'] > 0, 1, 0)

# --- read in riley_2011_30m_4326 = catalog.get_dataset('2011-climate-run-30m-4326').to_xarray()[['BP']] # there are slight mismatches with the coordinates in `ds` # interpolate riley data to exactly match ds coordinates riley_interp = riley_2011_30m_4326.interp( latitude=ds.latitude, longitude=ds.longitude, method='nearest' ) # assign coordinates ds['riley_BP_2011'] = riley_interp.BP # create burnable mask ds['riley_burnable_mask'] = xr.where(ds['riley_BP_2011'] > 0, 1, 0)

Historical fire perimeters¶

+ Filter the data¶

The Inter Agency Fire Perimeter History dataset includes small fires (where reported fire size limits are set by each reporting agency) and prescribed burns.

We filter out small fires in an effort to omit prescribed and more manageable fires. Based on the National Interagency Fire Center data, mean prescribed burns are around 30-50 acres (total fires / total acres). To omit most of these burns, we filter out all fire perimeters with an area less than 75 acres. (Note that there does not appear to be an input marking prescribed burns in the data).

In [4]:

# --- read in historical fire perimeter data (note, it can take a couple mins to load)
fp_path = 's3://carbonplan-ocr/evaluation/'
fp_name = 'InterAgencyFirePerimeterHistory_All_Years_View_-104997095188071827.gpkg'
gdf = gpd.read_file(os.path.join(fp_path, fp_name))
# convert CRS
gdf = gdf.to_crs('EPSG:4326')

# --- read in historical fire perimeter data (note, it can take a couple mins to load) fp_path = 's3://carbonplan-ocr/evaluation/' fp_name = 'InterAgencyFirePerimeterHistory_All_Years_View_-104997095188071827.gpkg' gdf = gpd.read_file(os.path.join(fp_path, fp_name)) # convert CRS gdf = gdf.to_crs('EPSG:4326')

In [5]:

# --- filter data
min_acres = 75
gdf = gdf[gdf['GIS_ACRES'] > min_acres]

# --- filter data min_acres = 75 gdf = gdf[gdf['GIS_ACRES'] > min_acres]

Add burn sum and mask and merge with ds¶

In [6]:

# --- define transform for rasterization
# Assuming lon/lat are 1D arrays
transform = rasterio.transform.from_bounds(
    ds.longitude.min(),
    ds.latitude.min(),
    ds.longitude.max(),
    ds.latitude.max(),
    len(ds.longitude),
    len(ds.latitude),
)

# --- compute burn sum
# create (geometry, value) tuples
shapes = [(geom, 1) for geom in gdf.geometry]

# burn all at once
burn_sum = rasterio.features.rasterize(
    [(geom, 1) for geom in gdf.geometry],
    out_shape=(len(ds.latitude), len(ds.longitude)),
    transform=transform,
    all_touched=True,
    dtype=np.uint16,  # allows counts > 255
    merge_alg=rasterio.features.MergeAlg.add,  # sum overlapping polygons
)

# add to dataset
ds['burn_sum'] = (('latitude', 'longitude'), burn_sum)

# get burn mask
ds['burn_mask'] = xr.where(ds['burn_sum'] > 0, 1, 0)

# --- define transform for rasterization # Assuming lon/lat are 1D arrays transform = rasterio.transform.from_bounds( ds.longitude.min(), ds.latitude.min(), ds.longitude.max(), ds.latitude.max(), len(ds.longitude), len(ds.latitude), ) # --- compute burn sum # create (geometry, value) tuples shapes = [(geom, 1) for geom in gdf.geometry] # burn all at once burn_sum = rasterio.features.rasterize( [(geom, 1) for geom in gdf.geometry], out_shape=(len(ds.latitude), len(ds.longitude)), transform=transform, all_touched=True, dtype=np.uint16, # allows counts > 255 merge_alg=rasterio.features.MergeAlg.add, # sum overlapping polygons ) # add to dataset ds['burn_sum'] = (('latitude', 'longitude'), burn_sum) # get burn mask ds['burn_mask'] = xr.where(ds['burn_sum'] > 0, 1, 0)

Save result¶

In [7]:

# --- delete everything we don't need
del burn_sum, shapes, transform, gdf, riley_interp, riley_2011_30m_4326
gc.collect()

# --- delete everything we don't need del burn_sum, shapes, transform, gdf, riley_interp, riley_2011_30m_4326 gc.collect()

Out[7]:

0

In [8]:

# --- rechunk
for name, var in ds.data_vars.items():
    ds[name] = var.chunk({'latitude': 6000, 'longitude': 4500})

# --- rechunk for name, var in ds.data_vars.items(): ds[name] = var.chunk({'latitude': 6000, 'longitude': 4500})

In [9]:

s3_path = 's3://carbonplan-ocr/evaluation/'
savename = 'benchmarking-input-dat.zarr'

# write to S3
ds.to_zarr(
    os.path.join(s3_path, savename),
    mode='w',  # overwrite (use "a" to append)
    compute=True,
    storage_options={'anon': False},  # set to True if public bucket
)

s3_path = 's3://carbonplan-ocr/evaluation/' savename = 'benchmarking-input-dat.zarr' # write to S3 ds.to_zarr( os.path.join(s3_path, savename), mode='w', # overwrite (use "a" to append) compute=True, storage_options={'anon': False}, # set to True if public bucket )

/opt/coiled/env/lib/python3.13/site-packages/zarr/api/asynchronous.py:244: ZarrUserWarning: Consolidated metadata is currently not part in the Zarr format 3 specification. It may not be supported by other zarr implementations and may change in the future.
  warnings.warn(

Out[9]:

<xarray.backends.zarr.ZarrStore at 0xf0c17606fa60>

In [ ]:

# ------------------------

# ------------------------

Compute quartiles for proportional plot¶

Restart kernel

Note that we omit burn probabilities of zero for this analysis because east of the 98th meridian BP=0 accounts for more than 20% of all data.

In [1]:

import gc
import os

import geopandas as gpd
import icechunk
import numpy as np
import rasterio.features

# import seaborn as sns    # plotting
import s3fs
import xarray as xr

from ocr import catalog  # for the riley data

# --- read in dat
s3_path = 's3://carbonplan-ocr/evaluation/'
savename = 'benchmarking-input-dat.zarr'
ds = xr.open_zarr(os.path.join(s3_path, savename))

import gc import os import geopandas as gpd import icechunk import numpy as np import rasterio.features # import seaborn as sns # plotting import s3fs import xarray as xr from ocr import catalog # for the riley data # --- read in dat s3_path = 's3://carbonplan-ocr/evaluation/' savename = 'benchmarking-input-dat.zarr' ds = xr.open_zarr(os.path.join(s3_path, savename))

Unclosed client session
client_session: <aiohttp.client.ClientSession object at 0xf534fff0ead0>
Unclosed connector
connections: ['deque([(<aiohttp.client_proto.ResponseHandler object at 0xf534ffe26390>, 151.370221169)])']
connector: <aiohttp.connector.TCPConnector object at 0xf534fff0ca50>
Unclosed client session
client_session: <aiohttp.client.ClientSession object at 0xf534fff0c7d0>
Unclosed connector
connections: ['deque([(<aiohttp.client_proto.ResponseHandler object at 0xf534ffe265d0>, 151.372040247)])']
connector: <aiohttp.connector.TCPConnector object at 0xf534fff0c410>
Unclosed client session
client_session: <aiohttp.client.ClientSession object at 0xf534fff0ed50>
Unclosed connector
connections: ['deque([(<aiohttp.client_proto.ResponseHandler object at 0xf534ffe26810>, 151.376278934)])']
connector: <aiohttp.connector.TCPConnector object at 0xf536b41c3390>

In [2]:

# --- filter to four slices
#     (CONUS, west, east, testbox) and later (all data, non-burnable)
slicenames = {
    'testbox': {'minlat': 42.7, 'maxlat': 46.3, 'minlon': -116.8, 'maxlon': -112.8},
    'CONUS': None,
    'West of -98': {'maxlon': -98},
    'East of -98': {'minlon': -98},
}

outdict = {}  # to hold outputs
for slc_key, slc_val in slicenames.items():
    if slc_key == 'testbox':
        outdict[slc_key] = ds.sel(
            latitude=slice(slc_val['maxlat'], slc_val['minlat']),
            longitude=slice(slc_val['minlon'], slc_val['maxlon']),
        )
    elif slc_key == 'CONUS':
        outdict[slc_key] = ds.copy()

    elif slc_key == 'West of -98':
        outdict[slc_key] = ds.where(ds['longitude'] < slc_val['maxlon'], drop=True)

    elif slc_key == 'East of -98':
        outdict[slc_key] = ds.where(ds['longitude'] >= slc_val['minlon'], drop=True)

# --- filter to four slices # (CONUS, west, east, testbox) and later (all data, non-burnable) slicenames = { 'testbox': {'minlat': 42.7, 'maxlat': 46.3, 'minlon': -116.8, 'maxlon': -112.8}, 'CONUS': None, 'West of -98': {'maxlon': -98}, 'East of -98': {'minlon': -98}, } outdict = {} # to hold outputs for slc_key, slc_val in slicenames.items(): if slc_key == 'testbox': outdict[slc_key] = ds.sel( latitude=slice(slc_val['maxlat'], slc_val['minlat']), longitude=slice(slc_val['minlon'], slc_val['maxlon']), ) elif slc_key == 'CONUS': outdict[slc_key] = ds.copy() elif slc_key == 'West of -98': outdict[slc_key] = ds.where(ds['longitude'] < slc_val['maxlon'], drop=True) elif slc_key == 'East of -98': outdict[slc_key] = ds.where(ds['longitude'] >= slc_val['minlon'], drop=True)

In [57]:

# --- FAIRLY MEMORY INTENSIVE (~400 GB...) ---
# --- split the dataset into 5 equal-area burn probability classes and compute expected vs observed burn proportions
#     assume each pixel has the same area
#
#     Note, this builds on the dataset created in Benchmark 1
#
save_on = True  # whether to save the data

for dsname, dsslice in outdict.items():
    # track progress
    print(f'Now solving: {dsname}')

    # pull out dat
    dsslice = outdict[dsname]
    bp = dsslice['burn_probability_2011']
    bp_nonzero = bp.where(bp > 0)
    burned = dsslice['burn_mask']

    # --- create quantile edges for 5 equal-area (equal-count) classes
    q_edges = np.linspace(0, 1, 6)  # 5 bins => 6 edges
    bp_quantiles = bp_nonzero.quantile(q_edges, dim=None).compute()
    edges = bp_quantiles.values

    # assign each pixel to a burn probability class using digitize
    # we broadcast bp_quantiles as edges
    def digitize_chunk(chunk, edges):
        return xr.DataArray(
            np.digitize(chunk, edges, right=True) - 1,
            dims=chunk.dims,
            coords=chunk.coords,
        )

    bp_class = xr.map_blocks(digitize_chunk, bp_nonzero, kwargs={'edges': edges})
    bp_class.name = 'bp_class'

    # --- bring bpclass into memory and get expected vs observed proportions
    total_bp_sum = bp_nonzero.sum()
    total_burned_sum = (burned == 1).sum()

    def masked_sum(var, cond):
        return var.where(cond).sum()

    expected_props = []
    observed_props = []
    # (for troubleshoot)
    expected_sum = []  # this should be the same value for all cases (observed sum can vary)

    for i in range(5):
        cond = (bp_nonzero >= edges[i]) & (bp_nonzero < edges[i + 1])
        expected_sum.append(cond.sum())
        expected_props.append(masked_sum(bp_nonzero, cond) / total_bp_sum)
        observed_props.append(masked_sum(burned == 1, cond) / total_burned_sum)

    expected_props = xr.concat(expected_props, dim='bp_class').compute()
    observed_props = xr.concat(observed_props, dim='bp_class').compute()
    expected_sum = xr.concat(expected_sum, dim='bp_class').compute()

    # check that expected bins are of equal size
    assert np.isclose(expected_sum.max(), expected_sum.min()), 'Expected sums are not equivalent'
    # check that both sum to 1
    assert np.isclose(expected_props.sum(), 1.0, atol=1e-02), 'Expected proportions do not sum to 1'
    assert np.isclose(observed_props.sum(), 1.0, atol=1e-02), 'Observed proportions do not sum to 1'

    # rename data and combine
    expected_props = expected_props.rename('bp2011_expected_burn_proportion')
    observed_props = observed_props.rename('bp2011_observed_burn_proportion')
    espected_sum = expected_sum.rename('bp2011_bpclass_sum')
    tmpds_props = xr.merge([expected_props, observed_props, expected_sum])
    tmpdf_props = tmpds_props.to_dataframe()
    tmpdf_props = tmpdf_props.reset_index()
    tmpdf_props['bp_class'] = tmpdf_props['bp_class'] + 1

    # save
    if save_on:
        s3 = s3fs.S3FileSystem()
        path = f's3://carbonplan-ocr/evaluation/benchmarking-processed/{dsname}_burnclass_proportions.parquet'
        with s3.open(path, 'wb') as f:
            tmpdf_props.to_parquet(f)

# --- FAIRLY MEMORY INTENSIVE (~400 GB...) --- # --- split the dataset into 5 equal-area burn probability classes and compute expected vs observed burn proportions # assume each pixel has the same area # # Note, this builds on the dataset created in Benchmark 1 # save_on = True # whether to save the data for dsname, dsslice in outdict.items(): # track progress print(f'Now solving: {dsname}') # pull out dat dsslice = outdict[dsname] bp = dsslice['burn_probability_2011'] bp_nonzero = bp.where(bp > 0) burned = dsslice['burn_mask'] # --- create quantile edges for 5 equal-area (equal-count) classes q_edges = np.linspace(0, 1, 6) # 5 bins => 6 edges bp_quantiles = bp_nonzero.quantile(q_edges, dim=None).compute() edges = bp_quantiles.values # assign each pixel to a burn probability class using digitize # we broadcast bp_quantiles as edges def digitize_chunk(chunk, edges): return xr.DataArray( np.digitize(chunk, edges, right=True) - 1, dims=chunk.dims, coords=chunk.coords, ) bp_class = xr.map_blocks(digitize_chunk, bp_nonzero, kwargs={'edges': edges}) bp_class.name = 'bp_class' # --- bring bpclass into memory and get expected vs observed proportions total_bp_sum = bp_nonzero.sum() total_burned_sum = (burned == 1).sum() def masked_sum(var, cond): return var.where(cond).sum() expected_props = [] observed_props = [] # (for troubleshoot) expected_sum = [] # this should be the same value for all cases (observed sum can vary) for i in range(5): cond = (bp_nonzero >= edges[i]) & (bp_nonzero < edges[i + 1]) expected_sum.append(cond.sum()) expected_props.append(masked_sum(bp_nonzero, cond) / total_bp_sum) observed_props.append(masked_sum(burned == 1, cond) / total_burned_sum) expected_props = xr.concat(expected_props, dim='bp_class').compute() observed_props = xr.concat(observed_props, dim='bp_class').compute() expected_sum = xr.concat(expected_sum, dim='bp_class').compute() # check that expected bins are of equal size assert np.isclose(expected_sum.max(), expected_sum.min()), 'Expected sums are not equivalent' # check that both sum to 1 assert np.isclose(expected_props.sum(), 1.0, atol=1e-02), 'Expected proportions do not sum to 1' assert np.isclose(observed_props.sum(), 1.0, atol=1e-02), 'Observed proportions do not sum to 1' # rename data and combine expected_props = expected_props.rename('bp2011_expected_burn_proportion') observed_props = observed_props.rename('bp2011_observed_burn_proportion') espected_sum = expected_sum.rename('bp2011_bpclass_sum') tmpds_props = xr.merge([expected_props, observed_props, expected_sum]) tmpdf_props = tmpds_props.to_dataframe() tmpdf_props = tmpdf_props.reset_index() tmpdf_props['bp_class'] = tmpdf_props['bp_class'] + 1 # save if save_on: s3 = s3fs.S3FileSystem() path = f's3://carbonplan-ocr/evaluation/benchmarking-processed/{dsname}_burnclass_proportions.parquet' with s3.open(path, 'wb') as f: tmpdf_props.to_parquet(f)

Now solving: testbox
Now solving: CONUS
Now solving: West of -98
Now solving: East of -98

In [3]:

# --- FAIRLY MEMORY INTENSIVE (~400 GB...) ---
#
# --- [ REPEAT FOR NON-BURNABLE ] ------------
#
# --- split the dataset into 5 equal-area burn probability classes and compute expected vs observed burn proportions
#     assume each pixel has the same area
#
#     Note, this builds on the dataset created in Benchmark 1
#
save_on = True  # whether to save the data

for dsname, dsslice in outdict.items():
    # track progress
    print(f'Now solving: {dsname}')

    # pull out dat
    dsslice = outdict[dsname]
    bp = dsslice['burn_probability_2011'].where(dsslice['riley_burnable_mask'] == 0)
    bp_nonzero = bp.where(bp > 0)
    burned = dsslice['burn_mask'].where(dsslice['riley_burnable_mask'] == 0)

    # --- create quantile edges for 5 equal-area (equal-count) classes
    q_edges = np.linspace(0, 1, 6)  # 5 bins => 6 edges
    bp_quantiles = bp_nonzero.quantile(q_edges, dim=None).compute()
    edges = bp_quantiles.values

    # assign each pixel to a burn probability class using digitize
    # we broadcast bp_quantiles as edges
    def digitize_chunk(chunk, edges):
        return xr.DataArray(
            np.digitize(chunk, edges, right=True) - 1,
            dims=chunk.dims,
            coords=chunk.coords,
        )

    bp_class = xr.map_blocks(digitize_chunk, bp_nonzero, kwargs={'edges': edges})
    bp_class.name = 'bp_class'

    # --- bring bpclass into memory and get expected vs observed proportions
    total_bp_sum = bp_nonzero.sum()
    total_burned_sum = (burned == 1).sum()

    def masked_sum(var, cond):
        return var.where(cond).sum()

    expected_props = []
    observed_props = []
    # (for troubleshoot)
    expected_sum = []  # this should be the same value for all cases (observed sum can vary)

    for i in range(5):
        cond = (bp_nonzero >= edges[i]) & (bp_nonzero < edges[i + 1])
        expected_sum.append(cond.sum())
        expected_props.append(masked_sum(bp_nonzero, cond) / total_bp_sum)
        observed_props.append(masked_sum(burned == 1, cond) / total_burned_sum)

    expected_props = xr.concat(expected_props, dim='bp_class').compute()
    observed_props = xr.concat(observed_props, dim='bp_class').compute()
    expected_sum = xr.concat(expected_sum, dim='bp_class').compute()

    # check that expected bins are of equal size
    assert np.isclose(expected_sum.max(), expected_sum.min()), 'Expected sums are not equivalent'
    # check that both sum to 1
    assert np.isclose(expected_props.sum(), 1.0, atol=1e-01), 'Expected proportions do not sum to 1'
    assert np.isclose(observed_props.sum(), 1.0, atol=1e-01), 'Observed proportions do not sum to 1'

    # rename data and combine
    expected_props = expected_props.rename('bp2011_expected_burn_proportion')
    observed_props = observed_props.rename('bp2011_observed_burn_proportion')
    espected_sum = expected_sum.rename('bp2011_bpclass_sum')
    tmpds_props = xr.merge([expected_props, observed_props, expected_sum])
    tmpdf_props = tmpds_props.to_dataframe()
    tmpdf_props = tmpdf_props.reset_index()
    tmpdf_props['bp_class'] = tmpdf_props['bp_class'] + 1

    # save
    if save_on:
        s3 = s3fs.S3FileSystem()
        path = f's3://carbonplan-ocr/evaluation/benchmarking-processed/{dsname}_burnclass_proportions_NB.parquet'
        with s3.open(path, 'wb') as f:
            tmpdf_props.to_parquet(f)

# --- FAIRLY MEMORY INTENSIVE (~400 GB...) --- # # --- [ REPEAT FOR NON-BURNABLE ] ------------ # # --- split the dataset into 5 equal-area burn probability classes and compute expected vs observed burn proportions # assume each pixel has the same area # # Note, this builds on the dataset created in Benchmark 1 # save_on = True # whether to save the data for dsname, dsslice in outdict.items(): # track progress print(f'Now solving: {dsname}') # pull out dat dsslice = outdict[dsname] bp = dsslice['burn_probability_2011'].where(dsslice['riley_burnable_mask'] == 0) bp_nonzero = bp.where(bp > 0) burned = dsslice['burn_mask'].where(dsslice['riley_burnable_mask'] == 0) # --- create quantile edges for 5 equal-area (equal-count) classes q_edges = np.linspace(0, 1, 6) # 5 bins => 6 edges bp_quantiles = bp_nonzero.quantile(q_edges, dim=None).compute() edges = bp_quantiles.values # assign each pixel to a burn probability class using digitize # we broadcast bp_quantiles as edges def digitize_chunk(chunk, edges): return xr.DataArray( np.digitize(chunk, edges, right=True) - 1, dims=chunk.dims, coords=chunk.coords, ) bp_class = xr.map_blocks(digitize_chunk, bp_nonzero, kwargs={'edges': edges}) bp_class.name = 'bp_class' # --- bring bpclass into memory and get expected vs observed proportions total_bp_sum = bp_nonzero.sum() total_burned_sum = (burned == 1).sum() def masked_sum(var, cond): return var.where(cond).sum() expected_props = [] observed_props = [] # (for troubleshoot) expected_sum = [] # this should be the same value for all cases (observed sum can vary) for i in range(5): cond = (bp_nonzero >= edges[i]) & (bp_nonzero < edges[i + 1]) expected_sum.append(cond.sum()) expected_props.append(masked_sum(bp_nonzero, cond) / total_bp_sum) observed_props.append(masked_sum(burned == 1, cond) / total_burned_sum) expected_props = xr.concat(expected_props, dim='bp_class').compute() observed_props = xr.concat(observed_props, dim='bp_class').compute() expected_sum = xr.concat(expected_sum, dim='bp_class').compute() # check that expected bins are of equal size assert np.isclose(expected_sum.max(), expected_sum.min()), 'Expected sums are not equivalent' # check that both sum to 1 assert np.isclose(expected_props.sum(), 1.0, atol=1e-01), 'Expected proportions do not sum to 1' assert np.isclose(observed_props.sum(), 1.0, atol=1e-01), 'Observed proportions do not sum to 1' # rename data and combine expected_props = expected_props.rename('bp2011_expected_burn_proportion') observed_props = observed_props.rename('bp2011_observed_burn_proportion') espected_sum = expected_sum.rename('bp2011_bpclass_sum') tmpds_props = xr.merge([expected_props, observed_props, expected_sum]) tmpdf_props = tmpds_props.to_dataframe() tmpdf_props = tmpdf_props.reset_index() tmpdf_props['bp_class'] = tmpdf_props['bp_class'] + 1 # save if save_on: s3 = s3fs.S3FileSystem() path = f's3://carbonplan-ocr/evaluation/benchmarking-processed/{dsname}_burnclass_proportions_NB.parquet' with s3.open(path, 'wb') as f: tmpdf_props.to_parquet(f)

Now solving: testbox
Now solving: CONUS
Now solving: West of -98
Now solving: East of -98

WIP: Compute LSS ...¶

Note: The random reference model is expensive and takes a long time to generate for CONUS (I haven't been patient enough to generate it yet). Consider just using the uniform model.

Compute the logarithmic score (LS) and logarithmic skill score (LSS)¶

Following Maron et al., 2025, the logarithmic score is defined by:

$\text{LS} = \frac{1}{N}\sum^{N}_{i-1} - \text{ln}\left(1-\left|p_i-o_i\right|\right)$

where $p_i$ is the burn probability of pixel $i$ and $o_i$ is the burn mask binary outcome (1 for burned, 0 for not burned).

We compute the $\text{LS}$ for the burn probability model as well as the two reference cases (random and uniform). THen we can compute the logarithmic skill score of the model relative to the reference cases by:

$\text{LSS} = 1 - \frac{\text{LS}}{\text{LS}_\text{ref}}$

where $\text{LS}_\text{ref}$ refers to either the random or uniform case logarithmic score.

In [3]:

import dask.array as da

dsslice = outdict['CONUS']


bp = dsslice['burn_probability_2011']
burn_mask = dsslice['burn_mask']

# --- reference models
# print("working on RANDOM model ... ")
# flat = bp.data.compute().flatten()           # compute into RAM
# np.random.shuffle(flat)                       # in-place shuffle
# shuffled = flat.reshape(bp.shape)
# bp_ref_random = xr.DataArray(shuffled, dims=bp.dims, coords=bp.coords, name="burn_probability_2011_REF_random")
# dsslice['burn_probability_2011_REF_random'] = bp_ref_random

# uniform mean
print('working on UNIFORM model ... ')
mean_bp = bp.mean()  # lazy if Dask
bp_ref_uniform = xr.full_like(bp, fill_value=mean_bp)
bp_ref_uniform.name = 'burn_probability_2011_REF_uniform'
dsslice['burn_probability_2011_REF_uniform'] = bp_ref_uniform

import dask.array as da dsslice = outdict['CONUS'] bp = dsslice['burn_probability_2011'] burn_mask = dsslice['burn_mask'] # --- reference models # print("working on RANDOM model ... ") # flat = bp.data.compute().flatten() # compute into RAM # np.random.shuffle(flat) # in-place shuffle # shuffled = flat.reshape(bp.shape) # bp_ref_random = xr.DataArray(shuffled, dims=bp.dims, coords=bp.coords, name="burn_probability_2011_REF_random") # dsslice['burn_probability_2011_REF_random'] = bp_ref_random # uniform mean print('working on UNIFORM model ... ') mean_bp = bp.mean() # lazy if Dask bp_ref_uniform = xr.full_like(bp, fill_value=mean_bp) bp_ref_uniform.name = 'burn_probability_2011_REF_uniform' dsslice['burn_probability_2011_REF_uniform'] = bp_ref_uniform

working on UNIFORM model ... 

In [4]:

dsslice['burn_probability_2011_REF_uniform']

dsslice['burn_probability_2011_REF_uniform']

Out[4]:

<xarray.DataArray 'burn_probability_2011_REF_uniform' (latitude: 97579,
                                                       longitude: 208881)> Size: 82GB
dask.array<full_like, shape=(97579, 208881), dtype=float32, chunksize=(6000, 4500), chunktype=numpy.ndarray>
Coordinates:
  * latitude   (latitude) float64 781kB 52.48 52.48 52.48 ... 22.43 22.43 22.43
  * longitude  (longitude) float64 2MB -128.4 -128.4 -128.4 ... -64.05 -64.05

xarray.DataArray

'burn_probability_2011_REF_uniform'

• latitude: 97579
• longitude: 208881

• dask.array<chunksize=(6000, 4500), meta=np.ndarray>

  Array Chunk Bytes 75.93 GiB 103.00 MiB Shape (97579, 208881) (6000, 4500) Dask graph 799 chunks in 1 graph layer Data type float32 numpy.ndarray

• Coordinates: (2)
  • latitude
    (latitude)
    float64
    52.48 52.48 52.48 ... 22.43 22.43

    units :

    degrees_north

    resolution :

    -0.0003079938784362739

    crs :

    EPSG:4326

    array([52.481695, 52.481387, 52.481079, ..., 22.428884, 22.428576, 22.428268],
          shape=(97579,))

  • longitude
    (longitude)
    float64
    -128.4 -128.4 ... -64.05 -64.05

    units :

    degrees_east

    resolution :

    0.0003079938784362739

    crs :

    EPSG:4326

    array([-128.387402, -128.387094, -128.386786, ...,  -64.054257,  -64.053949,
            -64.053641], shape=(208881,))

In [ ]:

# ... maybe this approach is faster...?


dsslice = outdict['CONUS']

bp = dsslice['burn_probability_2011']
burn_mask = dsslice['burn_mask']


bp_da = bp.data if isinstance(bp.data, da.Array) else da.from_array(bp.data, chunks=bp.chunks)
shape = bp_da.shape
# flatten lazily
flat = bp_da.reshape(-1)
# generate a permutation of indices lazily
idx = da.arange(flat.size)
idx = da.random.permutation(idx)  # note: this works lazily without chunks argument
# index into the flattened array lazily
shuffled_flat = flat[idx]
# reshape back to 2D
shuffled_bp_da = shuffled_flat.reshape(shape)
# wrap back in DataArray
bp_ref_random = xr.DataArray(
    shuffled_bp_da, dims=bp.dims, coords=bp.coords, name='burn_probability_2011_REF_random'
)

# ... maybe this approach is faster...? dsslice = outdict['CONUS'] bp = dsslice['burn_probability_2011'] burn_mask = dsslice['burn_mask'] bp_da = bp.data if isinstance(bp.data, da.Array) else da.from_array(bp.data, chunks=bp.chunks) shape = bp_da.shape # flatten lazily flat = bp_da.reshape(-1) # generate a permutation of indices lazily idx = da.arange(flat.size) idx = da.random.permutation(idx) # note: this works lazily without chunks argument # index into the flattened array lazily shuffled_flat = flat[idx] # reshape back to 2D shuffled_bp_da = shuffled_flat.reshape(shape) # wrap back in DataArray bp_ref_random = xr.DataArray( shuffled_bp_da, dims=bp.dims, coords=bp.coords, name='burn_probability_2011_REF_random' )

In [10]:

# dsslice['burn_probability_2011'].where('riley_burnable_mask' == 0)
dsslice['riley_burnable_mask'].where('riley_burnable_mask' == 0).plot()

# dsslice['burn_probability_2011'].where('riley_burnable_mask' == 0) dsslice['riley_burnable_mask'].where('riley_burnable_mask' == 0).plot()

Out[10]:

<matplotlib.collections.QuadMesh at 0xf9f2ccdbccd0>

In [13]:

observed_props.sum()

observed_props.sum()

Out[13]:

<xarray.DataArray 'burn_mask' ()> Size: 8B
array(0.95375784)

xarray.DataArray

'burn_mask'

• 0.9538

  array(0.95375784)

In [ ]:

In [ ]: