InsideDarkWeb.com

Large-scale cloud masking, grouping and mosaic in Descartes Labs platform

I am looking to create a stack of cloud-free sentinel-2 imagery by day for a large portion of California for all of 2020 in the Descartes Lab platform. I realize that this is a very large amount of imagery and will require significant memory to store and so I’m looking for methods to mosaic imagery by day. I’m new to Descartes and I think using the groupby and mosaic methods could work but not sure how to implement the groupby function on a masked image collection.

First, I establish a small area of interest in Northern California and dates from January to March of 2020 (I am using a small area and small time range to test out the workflow):

import descarteslabs as dl
import numpy as np

# Define a bounding box in N CA, in a GeoJSON
aoi = {
    "type": "Polygon",
    "coordinates": [
        [
            [-121.8430, 38.9676], 
            [-121.6991, 38.9676], 
            [-121.6991, 38.8590], 
            [-121.8430, 38.8590], 
        ]
    ],
}

start_datetime = "2019-01-01"
end_datetime = "2019-03-30"

Next, I create a Sentinel-2 scene stack with desired bands:

sentinel_scenes, sentinel_ctx = dl.scenes.search(
    aoi,
    products="sentinel-2:L1C",
    start_datetime = start_datetime,
    end_datetime = end_datetime,
#    cloud_fraction=0.7,
    limit = 500,
)

sentinel_all_stack = sentinel_scenes.stack('red green blue red-edge nir swir1 derived:ndvi', 
                                           sentinel_ctx, 
                                           processing_level='surface')

Then I create a Descartes Labs Sentinel-2 cloudmask stack for masking:

dlcloud_scenes, dlcloud_ctx = dl.scenes.search(
    aoi,
    products=["sentinel-2:L1C:dlcloud:v1"],
    start_datetime=start_datetime, end_datetime=end_datetime,
    limit=None
)


dlcloud_valid_cloudfree_stack = dlcloud_scenes.stack('valid_cloudfree', 
                                                     sentinel_ctx , 
                                                     data_type='Byte')

Next, I mask the Sentinel-2 collection with the Descartes Labs cloudmask:

dlcloud_valid_cloudfree_stack_rep = np.repeat(a=dlcloud_valid_cloudfree_stack,
                                          repeats=7,
                                          axis=1)

# Mask the pixel values where valid_cloud_free is 0
cloudfree_stack = np.ma.masked_where(dlcloud_valid_cloudfree_stack_rep==0, sentinel_all_stack)

From here, I would like to group the imagery by day. I tried:

for (year, month, day), day_scenes in imagery.groupby(
    "properties.date.year", "properties.date.month", "properties.date.day"):
    print("{}: {} scenes".format(month, np.array(day_scenes.stack("blue green red red-edge nir swir1 derived:ndvi", ctx)).shape))

And I get the error: ‘MaskedArray’ object has no attribute ‘groupby’

I would like to mosaic all images grouped by same day afterwards. All ideas are very welcome!

Geographic Information Systems Asked on November 14, 2021

1 Answers

One Answer

I believe the issue here may be that you are mixing the Workflows API (which has a groupby method) and the Scenes API (which does not). Full disclosure, I am a Descartes Labs employee.

Workflows is a bit different than Scenes in that is uses lazy proxy objects to define a series of processing steps, but does not actually do any of the data processing until the .compute(aoi) method is called for an aoi. More details are available here.

Here is an example showing how to create the image stack you are interested in using Workflows:

First we create a geocontext from the geojson of your aoi, this is used later when we run the workflow using .compute().

import descarteslabs as dl
import descarteslabs.workflows as wf
import numpy as np

# Define a bounding box in N CA, in a GeoJSON
aoi_geo = {
    "type": "Polygon",
    "coordinates": [
        [
            [-121.8430, 38.9676], 
            [-121.6991, 38.9676], 
            [-121.6991, 38.8590], 
            [-121.8430, 38.8590], 
        ]
    ],
}

aoi = wf.GeoContext(
    geometry=aoi_geo,
    resolution=10.0,
    crs='EPSG:3857',
    bounds_crs='EPSG:4326')

start_datetime = "2019-01-01"
end_datetime = "2019-02-28"

Here we define the Sentinel 2 imagery stack and cloud mask, then mask out cloudy pixels in each image.

s2_stack = (wf.ImageCollection.from_id('sentinel-2:L1C', 
                               start_datetime=start_datetime, 
                               end_datetime=end_datetime)
            .pick_bands('red green blue red-edge nir swir1 derived:ndvi')
           )

s2_cloud_mask = (wf.ImageCollection.from_id('sentinel-2:L1C:dlcloud:v1', 
                               start_datetime=start_datetime, 
                               end_datetime=end_datetime)
                .pick_bands('valid_cloudfree')
                )

s2_stack = s2_stack.concat_bands(s2_cloud_mask)

s2_masked = s2_stack.map(lambda img: img.mask(img.pick_bands('valid_cloudfree')==0))

Then we use the groupby and mosaic methods to create the daily image mosaic.

s2_daily = (s2_masked
            .groupby(dates=('year', 'month', 'day'))
            .mosaic()
            .pick_bands('red green blue red-edge nir swir1 derived:ndvi')
)

Finally we use the aoi to apply the data processing steps above to all the Sentinel 2 imagery in the aoi during the time period of interest.

s2_data = s2_daily.compute(aoi)

Answered by John Shriver on November 14, 2021

Add your own answers!

Related Questions

GeoServer Timestamp/Timezone error convert

2  Asked on September 30, 2021 by tim-van-heteren

       

Optimising a PostGIS table

2  Asked on September 29, 2021

     

Google Earth Engine – WFS, WMS, and KML feeds

1  Asked on September 29, 2021 by lane

       

Greyscale image exported

1  Asked on September 29, 2021

 

ST_AsPNG query has wrong position on Mapbox

0  Asked on September 29, 2021 by kawada

       

Export central coordinates of current map in print layout

1  Asked on September 29, 2021 by palemon

   

How to digitize building footprints with orthogonal edges?

3  Asked on September 28, 2021 by shawty

   

ArcMap cannot access WCS version 2.0.1 on MapServer?

1  Asked on September 28, 2021 by albert-tinon

       

Ask a Question

Get help from others!

© 2021 InsideDarkWeb.com. All rights reserved.