<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 be attributed to bothnatural factors (e.g. wildfire) and human activities (e.g. deforestation). Our approach is to buildupon an existing body of research that has provided high-resolution satellite-based estimates ofaboveground biomass (Spawn et al., 2020), forest cover change (Hansen et al., 2013), and changeattribution (Curtis et al., 2018). While many of the necessary data products already exist, we canintegrate, extend, or update these resources to provide global, current estimates that can beintegrated with the other resources produced by the coalition.Specifically, for any given spatial extent and time duration ($t1$ to $t2$), we can use threequantities — existing biomass, forest cover change, and change attribution — to estimate theeffective GHG emissions from land use changes. The simplest 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 spatial extent thatexperienced a stand-replacement disturbance between $t1$ and $t2$. The $TotalBiomass$ is estimatedas the aboveground biomass at time $t1$. This estimate can be further refined by attributing, foreach pixel, the source of forest cover loss (e.g. wildfire, deforestation, etc.), and using thosesources to express emissions fractionally and/or exclude certain categories from total estimates(e.g. rotational clear-cutting within tree plantations). Pixel-wise estimates can then be aggregatedinto 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 each direction (by summingemissions spatially), and then mosaic all the tiles together into a single xr dataset. Two notes forthis step. First, it gets a little funky if you're running using a client/cluster. So, it might beworth your time to just power down your cluster and run this on a single machine that you have forthis notebook. It'll run slower but your skull might be a little less bruised from banging your headagainst the wall. Second, you'll get a bunch of failed tiles below due to ocean cells that don'thave active tiles. That's great and normal. Bask in the glory of the the NaNs that will fill thatspace!<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 the year dimension.  Sanity checks! Let's take a look at the dataset and get a handle on whether this makes sense. Let's do it bylooking at a handful of big fires as a gut check for biomass loss. First let's check out the Thomas Fire in Santa Barbara in 2017, the second largest (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 this because there's anarrative here - we're on a Sherlock Holmes level hunt for fires. Using satellite biomass data asour 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 bigone, burning over 40 megatons of biomass in 2017. Some of that 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 further north. When we dothat, we notice that the 2017 Thomas Fire was actually dwarfed by something in 2007. What is thisaberration!? I thought wikipedia, the paragon of truth, told us that the Thomas Fire was thebiggest!<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 Santa Barbara county) in 2007.Over a million acres burned. And they were picked up by the loss in biomass. So this is somerejection of bias on my end! Excellent. Let's go check somewhere else just to see if it's working elsewhere. Let's travel up to my homestate of Washington State. As a gut check, the 2015 fire season was the largest in recent memory.Let's see if it jumps out at all in total 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 fire in 2015. So, let's seeif zooming in a little further into specific regions of big 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, one site of heavydeforestation 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 deforestation has 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 out Mexico since we canreasonably easily select a box around all of Mexico without including other countries. We includeGuatemala and a bit of Southern Texas/Arizona but let's just assume those are negligible compared tothe country of 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 of carbon!! :grin: our estimates are off by a factor of ~100. now, zarin's plot for mexico only includes deforestation,(so maybe that excludes the fire-attributed losses?). But they're def off by orders of magnitude.maybe these estimates we're calculating by combining everything is actually spatially-averagedquantities - 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 thinkthe coarsening effect might be okay as an approximation for now.HOWEVER, the per hectare unit thing is an issue- we summed things to go from 30m-->3km, so we'd needto divide them by ~10000 (to account for the coarsening 100-fold in each direction) to get back tothe per hectare unit. and then we use that per hectare unit and average it across regions andmultiply by the hectare area 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's area (in hectares) andthe units are on the same order of magnitude of zarin. now, this is very back-of-the-envelope. but ithink it makes sense! so, for reporting values, each of our coarsened pixels just needs to bedivided by 10,000 and then scaled by how many hectares that pixel is (which decreases going awayfrom 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 next episode, when wefigure out why our wiggles are totally off from theirs.... tune in!<jupyter_code>client.close()
cluster.close()
print("here!")<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11147857241284277, 'fraction_numerical': 0.029557668458040514, 'mean_word_length': 4.266690856313498, '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': '484474', 'license': 'MIT', 'n_tokens_mistral': 5160, 'n_tokens_neox': 4691, 'n_words': 1524, 'path': 'notebooks/v0/all_tiles.ipynb', 'repo_name': 'carbonplan/trace', 'size': '14526'}