Combining imagery, masking and sampling in Descartes Labs Platform

To run at scale, rather than calling .compute in a for-loop over each tile, you can submit all the tiles to the backend at once to run in parallel as Workflows Jobs, then asynchronously process them as they complete.

First though, you'll need to fix some edge cases with your cloud masking, and express your CDL crop masking in Workflows.

We'll walk through all of that below, but here's the same thing in notebook form, including interactive visualization and some other improvements which may be easier to follow.

Better cloud masking

l8_stack.concat_bands(l8_cloud_mask) assumes that l8_stack and l8_cloud_mask are the same length, which won't always be the case—some scenes might not have cloud masks processed. Here's a function to efficiently join the two, pick and mask only the scenes that have cloud masks, and construct the daily mosaics. It uses an ImageCollectionGroupby object (like a pandas groupby object) to efficiently lookup from the ImageCollections by date, and mosaic them at the same time:

def cloud_masked_daily_product(
    product_id: str, start_datetime: str, end_datetime: str
) -> wf.ImageCollection:
    "Get a product by ID, masked by the DL cloud mask and mosaicked by day"
    ic = wf.ImageCollection.from_id(product_id, start_datetime, end_datetime)
    cloudmask = (
        wf.ImageCollection.from_id(
            product_id + ":dlcloud:v1", start_datetime, end_datetime
        ).pick_bands("valid_cloudfree")
        == 0
    )

    # Make an ImageCollectionGroupby object, for quicker lookups
    # from `ic` by date (you can use it like a dict)
    ic_date_groupby = ic.groupby(dates=("year", "month", "day"))
    # For each cloudmask date, pick the corresponding image from `ic` by date, mosiac both, and mask them.
    # (Not all scenes have cloudmasks processed, so this ensures we only return scenes that do.)
    return cloudmask.groupby(dates=("year", "month", "day")).map(
        lambda ymd, mask_imgs: ic_date_groupby[ymd].mosaic().mask(mask_imgs.mosaic())
    )

We'll also write a quick NDVI function, rather than using the derived:ndvi band (which can be slower):

def ndvi(ic: wf.ImageCollection) -> wf.ImageCollection:
    nir, red = ic.unpack_bands("nir red")
    ndvi = (nir - red) / (nir + red)
    return ndvi.rename_bands("ndvi")

Putting these together, we can prepare our input imagery more concisely:

l8_daily = cloud_masked_daily_product(
    "landsat:LC08:01:T1:TOAR", start_datetime, end_datetime
).pick_bands("red green blue nir swir1")
l8_with_ndvi = l8_daily.concat_bands(ndvi(l8_daily))

s2_daily = cloud_masked_daily_product(
  "sentinel-2:L1C", start_datetime, end_datetime
).pick_bands("red green blue nir swir1")
s2_with_ndvi = s2_daily.concat_bands(ndvi(s2_daily))

s1 = wf.ImageCollection.from_id(
    "sentinel-1:GRD", start_datetime, end_datetime
).pick_bands("vh vv")
s1_daily = s1.groupby(dates=("year", "month", "day")).mosaic()

Incorporating CDL

Workflows doesn't have a built-in equivalent to np.isin, but we can write one. It's not quite as efficient, but the difference isn't noticeable.

def isin(ic: wf.ImageCollection, values: list) -> wf.ImageCollection:
    "Like np.isin, for Workflows"
    assert len(values) > 0
    result = False
    for value in values:
        result = result | (ic == value)
    return result

Using that, we'll replicate your CDL classification logic:

# picking all 4 years of CDL at once is more efficient
cdl = wf.ImageCollection.from_id(
    "usda:cdl:v1", start_datetime="2016-12-01", end_datetime="2020-01-01"
).pick_bands("class")

grains_oils_grass_beans = [1,2,3,4,5,6,10,11,12,13,21,22,23,24,25,26,27,28,29,
                        30,31,32,33,34,35,36,37,38,39,41,42,43,44,45,46,51,
                        52,53,225,226,228,230,232,234,235,236,237,238,239,240,241,254]

deli_crops = [14, 48, 49, 50, 54, 55, 57, 206, 207, 208, 209, 213, 214, 216,
            219, 221, 222, 224, 227, 229, 231, 242, 243, 244, 245, 246, 247,
            248, 249, 250]

tree_crops = [66, 67, 68, 69, 72, 74, 75, 76, 77, 204, 210, 211, 212, 215, 217,
            218,220, 223]

crops_list = deli_crops + tree_crops

is_crops = isin(cdl, crops_list)
is_crops_19 = is_crops[-1]

four_year_combo = is_crops.sum(axis="images") + is_crops_19  # double-weight 2019

# create binary array where 1 = cultivated for 2 years+ and/or in 2019,
# 0 = cultivated for only 1 year and not in 2019
four_year_binary = four_year_combo >= 2

# invert it to mask form, where True means masked
cdl_mask = ~four_year_binary

Then, we apply the cdl_mask to all of our training imagery:

l8_masked = l8_with_ndvi.mask(cdl_mask)
s2_masked = s2_with_ndvi.mask(cdl_mask)
s1_masked = s1_daily.mask(cdl_mask)

To fetch this data, you can compute all three ImageCollections at once, which is faster because it will happen in parallel, and the CDL mask will be reused for all three:

l8_data, s2_data, s1_data = wf.compute(
  [l8_masked.ndarray, s2_masked.ndarray, s1_masked.ndarray], tile
)

Now you have all of your input arrays, masked by CDL!

One thing to notice is that because the three satellites have different revisit rates, these three imagery stacks are different lengths. If your model requires having the same number of input scenes from each sensor, you'll need to add more logic to composite or filter the scenes down.

Also, I noticed that many of your tiles don't contain any of the crops you're looking for—they'd all be masked out. In the notebook I have code for filtering out these tiles ahead-of-time.

Running at scale

Submitting all the tiles at once will be much faster than calling .compute one-at-a-time, but also requires some boilerplate code.

We'll call .compute with block=False to immediately get a Job object, instead of waiting for the job to finish. Then, we'll write a function to check which Jobs are done, and process them in whatever order they finish.

# submit all the tiles to run, asynchronously
jobs = [
  wf.compute(
    [l8_masked.ndarray, s2_masked.ndarray, s1_masked.ndarray],
    tile,
    block=False
  )
  for tile in tqdm(tiles_to_run)
]

Then here's a generator function to iterate over the jobs as they complete:

from typing import Iterator, Sequence, Optional
import time

def as_completed(jobs: Sequence[wf.Job], interval_sec: Optional[int] = None) -> Iterator[wf.Job]:
    """
    Iterator over Jobs that yields each Job when it completes.
    
    Parameters
    ----------
    jobs: Sequence[wf.Job]
        Jobs to wait for
    interval_sec: Optional[int], default None
        Wait at least this many seconds between polling for job updates.
        If None (default), uses ``max(5, len(jobs) // 5)``.
        
    Yields
    ------
    job: wf.Job
        A completed job (either succeeded or failed).
    """
    jobs = list(jobs)
    if interval_sec is None:
        inverval_sec = max(5, len(jobs) // 5)
    while len(jobs) > 0:
        loop_start = time.perf_counter()

        i = 0
        while i < len(jobs):
            job = jobs[i]
            if not job.done:  # in case it's already loaded
                try:
                    job.refresh()
                except Exception:
                    continue  # be resilient to transient errors for now

            if job.done:
                yield job
                del jobs[i]  # "advances" i
            else:
                i += 1

        loop_duration = time.perf_counter() - loop_start
        if len(jobs) > 0 and loop_duration < interval_sec:
            time.sleep(interval_sec - loop_duration)

Finally, we'll use this function to iterate over the jobs list and handle each result:

failed = []
for job in as_completed(jobs):
    if job.error is not None:
        failed.append(job)
        print(job.error)
    else:
        l8_data, s2_data, s1_data = job.result(progress_bar=False)
        handle_result(l8_data, s2_data, s1_data)

Again, here is all of this in a working notebook, which also has examples of you'd visualize all the intermediate parts on an interactive map.