I am trying to sample data from Descartes Lab imagery to train a random forest classifier in the Descartes Plaform for an agricultural region. I have already written a script to prepare the imagery of interest (Sentinel-2, Sentinel-1 and Landsat 8) using the Descartes Labs cloud masks where applicable, and I have not figured out yet how to sample points. I see an example in the Descartes Labs platform that imports training data the example, but I would like to create my samples within the platform. How can I create a sample data set within the platform at scale?
First I import the packages of interest, identify my AOI (the central valley of California) and dates of interest:
# Import packages
import descarteslabs as dl
import descarteslabs.workflows as wf
import numpy as np
from concurrent.futures import ThreadPoolExecutor
from shapely.geometry import shape, box
from random import random
from tqdm.notebook import tqdm
sac = shape(dl.places.shape(
'north-america_united-states_california_sacramento-valley'
).geometry)
sj = shape(dl.places.shape(
'north-america_united-states_california_san-joaquin-valley'
).geometry)
central_valley_aoi = sac.union(sj)
start_datetime = "2019-01-01"
end_datetime = "2019-03-28"
I then create tiles with my AOI:
tiles = dl.raster.dltiles_from_shape(resolution=10,
tilesize=500,
pad=0,
shape=central_valley_aoi)
I then write out my functions to prep each Landsat 8, Sentinel-1 and Sentinel-2 imagery using . Landsat 8:
l8_stack = (wf.ImageCollection.from_id('landsat:LC08:01:T1:TOAR',
start_datetime=start_datetime,
end_datetime=end_datetime)
.pick_bands('red green blue nir swir1 derived:ndvi')
)
l8_cloud_mask = (wf.ImageCollection.from_id('landsat:LC08:01:T1:TOAR:dlcloud:v1',
start_datetime=start_datetime,
end_datetime=end_datetime)
.pick_bands('valid_cloudfree')
)
l8_stack = l8_stack.concat_bands(l8_cloud_mask)
l8_masked = l8_stack.map(lambda img: img.mask(img.pick_bands('valid_cloudfree')==0))
l8_daily = (l8_masked
.groupby(dates=('year', 'month', 'day'))
.mosaic()
.pick_bands('red green blue nir swir1 derived:ndvi')
)
Sentinel-2:
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))
s2_daily = (s2_masked
.groupby(dates=('year', 'month', 'day'))
.mosaic()
.pick_bands('red green blue red-edge nir swir1 derived:ndvi')
)
Sentinel-1:
s1_stack = (wf.ImageCollection.from_id('sentinel-1:GRD',
start_datetime=start_datetime,
end_datetime=end_datetime)
.pick_bands('vh vv')
)
s1_daily = (s1_stack
.groupby(dates=('year', 'month', 'day'))
.mosaic()
.pick_bands('vh vv' )
)
These functions to prep each satellite imagery can be run easily with these simple lines and ideally will be called by tile, although I haven't figured out how to integrate the Cropland Mask:
s2_data = s2_daily.compute(tile)
l8_data = l8_daily.compute(tile)
s1_data = s1_daily.compute(tile)
I then prep the cropland mask (where specific crops have been grown in California more than twice in the past for year and/or in 2019 to create a binary array:
cropland_2016, ctx = dl.scenes.search(tile_cord,
products='usda:cdl:v1',
start_datetime="2016-12-01",
end_datetime="2017-01-01",
limit=5
)
cropland_2017, ctx = dl.scenes.search(tile_cord,
products='usda:cdl:v1',
start_datetime="2017-12-01",
end_datetime="2018-01-01",
limit=5
)
cropland_2018, ctx = dl.scenes.search(tile_cord,
products='usda:cdl:v1',
start_datetime="2018-12-01",
end_datetime="2019-01-01",
limit=5
)
cropland_2019, ctx = dl.scenes.search(tile_cord,
products='usda:cdl:v1',
start_datetime="2019-12-01",
end_datetime="2020-01-01",
limit=5
)
#Arrays of Cropland Data Layer by year
cld_16 = cropland_2016[0].ndarray('class', ctx)
cld_17 = cropland_2017[0].ndarray('class', ctx)
cld_18 = cropland_2018[0].ndarray('class', ctx)
cld_19 = cropland_2019[0].ndarray('class', ctx)
#Cropland Data Layer Codes, by crop group
#https://www.nass.usda.gov/Research_and_Science/Cropland/metadata/metadata_ca19.htm
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
# binary remapping of Cropland Data Layer to include only Delicate Crops
cld_16_deli = np.isin(cld_16, [crops_list]).astype(int)
cld_17_deli = np.isin(cld_17, [crops_list]).astype(int)
cld_18_deli = np.isin(cld_18, [crops_list]).astype(int)
cld_19_deli = np.isin(cld_19, [crops_list]).astype(int)
# weighting 2019 double
cld_19_deli2 = cld_19_deli*2
# combine all years of the Cropland Data Layer binary array
four_year_combo = cld_16_deli + cld_17_deli + cld_18_deli + cld_19_deli2
# 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 = np.isin(four_year_combo, [2,3,4,5]).astype(int)
I would then like to call the imagery collection for each tile and mask all pixels using the corresponding cropland data layer mask. I have not finalized this piece.
# something like this?
tile_data = {}
for k, tile in tqdm(enumerate(tiles['features'])):
s2_data = s2_daily.compute(tile)
l8_data = l8_daily.compute(tile)
s1_data = s1_daily.compute(tile)
# stack all these bands
# mask with cropland data layer to keep only areas that have grown desired crops
}
At this point, once I have my stack of masked imagery, I would like to extract training data from the 2019 imagery using the 2019 Cropland Data Layer classification. How would I approach this in the Descartes Labs Platform at scale?