<jupyter_start><jupyter_text>Forest Emissions Tracking - Phase I_by <NAME> and <NAME> (CarbonPlan)_March 29, 2020 IntroductionIn general, greenhouse gasses (GHGs) arising from forest land use changes can beattributed to both natural factors (e.g. wildfire) and human activities (e.g.deforestation). Our approach is to build upon an existing body of research thathas provided high-resolution satellite-based estimates of aboveground biomass(Spawn et al., 2020), forest cover change (Hansen et al., 2013), and changeattribution (Curtis et al., 2018). While many of the necessary data productsalready exist, we can integrate, extend, or update these resources to provideglobal, current estimates that can be integrated with the other resourcesproduced by the coalition.Specifically, for any given spatial extent and time duration ($t1$ to $t2$), wecan use three quantities — existing biomass, forest cover change, and changeattribution — to estimate the effective GHG emissions from land use changes. Thesimplest estimate is:$\Delta Biomass (t) = TotalBiomass (t) * \Delta ForestCover (\%)$$Emissions (tCO_2) = \Delta Biomass (t) * 0.5 (tC/t) * 3.67 (tC02 / tC)$where $\Delta ForestCover$ is the fraction of pixels within the given spatialextent that experienced a stand-replacement disturbance between $t1$ and $t2$.The $TotalBiomass$ is estimated as the aboveground biomass at time $t1$. Thisestimate can be further refined by attributing, for each pixel, the source offorest cover loss (e.g. wildfire, deforestation, etc.), and using those sourcesto express emissions fractionally and/or exclude certain categories from totalestimates (e.g. rotational clear-cutting within tree plantations). Pixel-wiseestimates can then be aggregated into province and country-wide estimates.  SetupTo begin, we'll import a handful of Python libraries and set a few constants.<jupyter_code>%matplotlib inline
# import dask
import intake
import xarray as xr
from tqdm import tqdm
import numcodecs
import numpy as np

TC02_PER_TC = 3.67
TC_PER_TBM = 0.5
SQM_PER_HA = 10000
ORNL_SCALING = 0.1
from dask_gateway import Gateway

gateway = Gateway()
options = gateway.cluster_options()
options.worker_cores = 2
options.worker_memory = 24
cluster = gateway.new_cluster(cluster_options=options)
cluster.adapt(minimum=1, maximum=300)
cluster
client = cluster.get_client()
client
# data catalog
cat = intake.open_catalog(
    "https://raw.githubusercontent.com/carbonplan/forest-emissions-tracking/master/catalog.yaml"
)
import fsspec

with fsspec.open(
    "https://storage.googleapis.com/earthenginepartners-hansen/GFC-2018-v1.6/treecover2000.txt"
) as f:
    lines = f.read().decode().splitlines()
print(len(lines))

# americas tiles
lat_tags = [
    "80N",
    "70N",
    "60N",
    "50N",
    "40N",
    "30N",
    "20N",
    "10N",
    "00N",
    "10S",
    "20S",
    "30S",
    "40S",
    "50S",
]
lon_tags = [f"{n:03}W" for n in np.arange(10, 190, 10)] + [
    f"{n:03}E" for n in np.arange(0, 190, 10)
]

# all tiles
lats = []
lons = []
for line in lines:
    pieces = line.split("_")
    lat = pieces[-2]
    lon = pieces[-1].split(".")[0]

    if (lat in lat_tags) and (lon in lon_tags):
        lats.append(lat)
        lons.append(lon)
print(len(lats))
def _preprocess(da, lat=None, lon=None):
    da = da.rename({"x": "lon", "y": "lat"}).squeeze(drop=True)
    if lat is not None:
        da = da.assign_coords(lat=lat, lon=lon)
    return da


def open_hansen_2018_tile(lat, lon, emissions=False):
    ds = xr.Dataset()

    # Min Hansen data
    variables = [
        "treecover2000",
        "gain",
        "lossyear",
        "datamask",
    ]  # , "first", "last"]
    for v in variables:
        da = (
            cat.hansen_2018(variable=v, lat=lat, lon=lon)
            .to_dask()
            .pipe(_preprocess)
        )
        # force coords to be identical
        if ds:
            da = da.assign_coords(lat=ds.lat, lon=ds.lon)
        ds[v] = da

    ds["treecover2000"] /= 100.0
    ds["lossyear"] += 2000

    # Hansen biomass
    ds["agb"] = (
        cat.hansen_biomass(lat=lat, lon=lon)
        .to_dask()
        .pipe(_preprocess, lat=ds.lat, lon=ds.lon)
    )
    if emissions:
        # Hansen emissions
        ds["emissions_ha"] = (
            cat.hansen_emissions_ha(lat=lat, lon=lon)
            .to_dask()
            .pipe(_preprocess, lat=ds.lat, lon=ds.lon)
        )
        ds["emissions_px"] = (
            cat.hansen_emissions_px(lat=lat, lon=lon)
            .to_dask()
            .pipe(_preprocess, lat=ds.lat, lon=ds.lon)
        )

    return ds
# open a single 10x10degree tile of the Hansen 30x30m data
lat = lats[1]
lon = lons[1]
box = dict(lat=slice(0, 40000, 100), lon=slice(0, 40000, 100))

ds = open_hansen_2018_tile(lat, lon)
display(ds)
encoding = {"emissions": {"compressor": numcodecs.Blosc()}}
@dask.delayed
def process_one_tile(lat, lon):
    url = f"gs://carbonplan-scratch/global-forest-emissions/{lat}_{lon}.zarr"

    mapper = fsspec.get_mapper(url)

    with dask.config.set(scheduler="threads"):
        ds = open_hansen_2018_tile(lat, lon)
        ds = calc_one_tile(ds)[["emissions"]]
        ds = ds.chunk({"lat": 4000, "lon": 4000, "year": 2})
        ds.to_zarr(mapper, encoding=encoding, mode="w")
        return url
def compute_grid_area(da):
    R = 6.371e6
    total_area = 4 * np.pi * R ** 2
    dϕ = np.radians((da["lat"][1] - da["lat"][0]).values)
    dλ = np.radians((da["lon"][1] - da["lon"][0]).values)
    dA = R ** 2 * np.abs(dϕ * dλ) * np.cos(np.radians(da["lat"]))
    areacella = dA * (0 * da + 1)
    return areacella / SQM_PER_HA
from itertools import product

all_tiles = list(product(lat_tags, lon_tags))
rerun_tiles = False
if rerun_tiles:
    tiles = []
    for lat, lon in tqdm(list(zip(lats, lons))):
        tiles.append(client.persist(process_one_tile(lat, lon), retries=1))<jupyter_output><empty_output><jupyter_text>Next, we'll load each of the tiles from above, coarsen them 100-fold in eachdirection (by summing emissions spatially), and then mosaic all the tilestogether into a single xr dataset. Two notes for this step. First, it gets alittle funky if you're running using a client/cluster. So, it might be worthyour time to just power down your cluster and run this on a single machine thatyou have for this notebook. It'll run slower but your skull might be a littleless bruised from banging your head against the wall. Second, you'll get a bunchof failed tiles below due to ocean cells that don't have active tiles. That'sgreat and normal. Bask in the glory of the the NaNs that will fill that space!<jupyter_code>len(list(zip(lats, lons)))
recombine = True
if recombine:
    list_all_coarsened = []
    for lat, lon in tqdm(list(zip(lats, lons))):
        try:
            mapper = fsspec.get_mapper(
                f"gs://carbonplan-scratch/global-forest-emissions/{lat}_{lon}.zarr"
            )
            da_global = xr.open_zarr(mapper)
            da_mask = da_global.isel(year=0, drop=True)
            da_area = compute_grid_area(da_mask)
            list_all_coarsened.append(
                (da_global * da_area).coarsen(lat=100, lon=100).sum().compute()
            )
        except:
            print("{} {} didnt work booooo".format(lat, lon))<jupyter_output><empty_output><jupyter_text>Then combine that whole list into a single dataset and write it out.<jupyter_code>coarsened_url = (
    f"gs://carbonplan-scratch/global-forest-emissions/global/3000m/raster.zarr"
)
mapper = fsspec.get_mapper(coarsened_url)
recombine = False

if recombine:
    # with dask.config.set(**{'array.slicing.split_large_chunks': False}):
    combined_ds = xr.combine_by_coords(
        list_all_coarsened, compat="override", coords="minimal"
    )
#     combined_ds.to_zarr(mapper,
#                      encoding=encoding,
#                      mode='w')
else:
    combined_ds = xr.open_zarr(mapper)
#     combined_ds = combined_ds.rename({'emissions': 'Emissions [tC02/year]'})
combined_ds.coarsen(lat=100, lon=100).isel(year=0).plot()
coarsened_sel = combined_ds.isel(year=0).coarsen(lat=10, lon=10).sum().compute()
coarsened_sel.emissions.plot()
combined_ds.isel(year=17).emissions.plot()
diff = combined_ds_testing.sel(year=2005) - combined_ds.sel(year=2005).sel(
    lat=combined_ds_testing.lat, lon=combined_ds_testing.lon
)
diff.emissions.plot()
combined_ds = combined_ds.rename({"emissions": "Emissions [tC02/year]"})<jupyter_output><empty_output><jupyter_text>Then we'll load all of these back into a single array and concatenate along theyear dimension.  Sanity checks! Let's take a look at the dataset and get a handle on whether this makes sense.Let's do it by looking at a handful of big fires as a gut check for biomassloss. First let's check out the Thomas Fire in Santa Barbara in 2017, the secondlargest (outside of 2020) California history since we don't have 2020 data yet.This ranking is according to wikipedia -https://en.wikipedia.org/wiki/List_of_California_wildfires. I include thisbecause there's a narrative here - we're on a Sherlock Holmes level hunt forfires. Using satellite biomass data as our clues!<jupyter_code>santa_barbara = combined_ds.sel(lat=slice(34, 34.5), lon=slice(-119.75, -119))[
    "Emissions [tC02/year]"
]
santa_barbara.sel(year=2017).plot()<jupyter_output><empty_output><jupyter_text>But then let's see how that fire compared to years previous. It turns out that,indeed, it was a big one, burning over 40 megatons of biomass in 2017. Some ofthat fire carried over into 2018 too.<jupyter_code>santa_barbara.sum(dim=["lat", "lon"]).plot()<jupyter_output><empty_output><jupyter_text>Out of curiosity, we can take a look at a slightly bigger box, extending furthernorth. When we do that, we notice that the 2017 Thomas Fire was actually dwarfedby something in 2007. What is this aberration!? I thought wikipedia, the paragonof truth, told us that the Thomas Fire was the biggest!<jupyter_code>santa_barbara_bigger = combined_ds.sel(
    lat=slice(34, 35), lon=slice(-119.75, -119)
)["Emissions [tC02/year]"]
santa_barbara_bigger.sum(dim=["lat", "lon"]).plot()<jupyter_output><empty_output><jupyter_text>Turns out that actually there was _giant_ complex of fires (also in SantaBarbara county) in 2007. Over a million acres burned. And they were picked up bythe loss in biomass. So this is some rejection of bias on my end! Excellent. Let's go check somewhere else just to see if it's working elsewhere. Let'stravel up to my home state of Washington State. As a gut check, the 2015 fireseason was the largest in recent memory. Let's see if it jumps out at all intotal Washington losses.<jupyter_code>washington_state = combined_ds.sel(
    lat=slice(45.5, 49), lon=slice(-124.6, -117)
)["Emissions [tC02/year]"]
washington_state.sum(dim=["lat", "lon"]).plot()<jupyter_output><empty_output><jupyter_text>Well, let's take a look at the Okanaogan region where they had a giant firein 2015. So, let's see if zooming in a little further into specific regions ofbig fires leap out.<jupyter_code>washington_state.sel(lat=slice(48.4, 48.6), lon=slice(-119.8, -119.6)).sum(
    dim=["lat", "lon"]
).plot()<jupyter_output><empty_output><jupyter_text>Now let's check out Brazil. First we can look at the state of Mato Grosso, onesite of heavy deforestation for conversion to agriculure.<jupyter_code>combined_ds.sel(lat=slice(-14, -11), lon=slice(-55, -54)).sum(
    dim=["lat", "lon"]
)["Emissions [tC02/year]"].plot()<jupyter_output><empty_output><jupyter_text>Then we can switch to the state of Pará.<jupyter_code>combined_ds.sel(lat=slice(-9, -2), lon=slice(-54, -49)).sum(dim=["lat", "lon"])[
    "Emissions [tC02/year]"
].plot()<jupyter_output><empty_output><jupyter_text>Selecting a year we can see individual plots of land where mass deforestationhas occurred.<jupyter_code>combined_ds.sel(lat=slice(-9, -2), lon=slice(-54, -49)).sel(year=2017)[
    "Emissions [tC02/year]"
].plot(vmax=1e6)<jupyter_output><empty_output><jupyter_text>Gutcheck by trying to reproduce Zarin et al (2016) a.k.a. Reproducibility Rumble Now let's try to reproduce a figure from Zarin et al (2016). We'll try outMexico since we can reasonably easily select a box around all of Mexico withoutincluding other countries. We include Guatemala and a bit of SouthernTexas/Arizona but let's just assume those are negligible compared to the countryof Mexico. So, comparing this to the plot from<jupyter_code>mexico = combined_ds.sel(lat=slice(14, 31), lon=slice(-117, -87))
mexico.sum(dim=["lat", "lon"])["emissions"].plot()<jupyter_output><empty_output><jupyter_text>SAD FACE! We're off by a _ton_ <--- a pun! because we're talking about tons ofcarbon!! :grin: our estimates are off by a factor of ~100. now, zarin's plot for mexico onlyincludes deforestation, (so maybe that excludes the fire-attributed losses?).But they're def off by orders of magnitude. maybe these estimates we'recalculating by combining everything is actually spatially-averaged quantities -tCO2/hectare. which would mean that an average is really what we want. and thatcoarsening isn't quite right because the geographic transformation matters.nevertheless, i think the coarsening effect might be okay as an approximationfor now.HOWEVER, the per hectare unit thing is an issue- we summed things to go from30m-->3km, so we'd need to divide them by ~10000 (to account for the coarsening100-fold in each direction) to get back to the per hectare unit. and then we usethat per hectare unit and average it across regions and multiply by the hectarearea to get the right tCO2. as a gut check i applied that "divide by 10,000"correction that around my bounding box for mexico, multiplied it by mexico'sarea (in hectares) and the units are on the same order of magnitude of zarin.now, this is very back-of-the-envelope. but i think it makes sense! so, forreporting values, each of our coarsened pixels just needs to be divided by10,000 and then scaled by how many hectares that pixel is (which decreases goingaway from the equator, but we could make a raster of that).<jupyter_code>hectares_in_mexico = 197253494
correction_for_having_summed_during_coarsening = 100 * 100
(
    mexico.mean(dim=["lat", "lon"])
    * hectares_in_mexico
    / correction_for_having_summed_during_coarsening
)["Emissions [tC02/year]"].plot()<jupyter_output><empty_output><jupyter_text>ba da bing ba da boom! Now we're in the ballpark of Zarin. Don't miss the nextepisode, when we figure out why our wiggles are totally off from theirs.... tunein!<jupyter_code>client.close()
cluster.close()
print("here!")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11121014741172437, 'fraction_numerical': 0.029413781282139183, 'mean_word_length': 4.166843783209352, 'pattern_counts': {'":': 8, '<': 60, '<?xml version=': 0, '>': 60, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '178337', 'license': 'MIT', 'n_tokens_mistral': 5194, 'n_tokens_neox': 4721, 'n_words': 1513, 'path': 'notebooks/all_tiles.ipynb', 'repo_name': 'orianac/forest-emissions-tracking', 'size': '14597'}