<jupyter_start><jupyter_text>ARDC Training: Python NotebooksTask-A: Cloud-free Mosaics and K-means Clustering >  Import the Datacube Configuration<jupyter_code>import datacube
import utils.data_cube_utilities.data_access_api as dc_api  

from datacube.utils.aws import configure_s3_access
configure_s3_access(requester_pays=True)

api = dc_api.DataAccessApi()
dc = datacube.Datacube(app = 'ardc_task_a')
api.dc = dc<jupyter_output><empty_output><jupyter_text>> Browse the available Data Cubes<jupyter_code>list_of_products = dc.list_products()
netCDF_products = list_of_products[list_of_products['format'] == 'NetCDF']
netCDF_products<jupyter_output><empty_output><jupyter_text>> Pick a product  >Use the platform and product names from the previous block to select a Data Cube.<jupyter_code># Change the data platform and data cube here

platform = 'LANDSAT_7'
product = 'ls7_usgs_sr_scene'<jupyter_output><empty_output><jupyter_text>>  Display Latitude-Longitude and Time Bounds of the Data Cube<jupyter_code>from utils.data_cube_utilities.dc_time import _n64_to_datetime, dt_to_str

extents = api.get_full_dataset_extent(platform = platform, product = product, measurements=[])

latitude_extents = (min(extents['latitude'].values),max(extents['latitude'].values))
longitude_extents = (min(extents['longitude'].values),max(extents['longitude'].values))
time_extents = (min(extents['time'].values),max(extents['time'].values))

print("Latitude Extents:", latitude_extents)
print("Longitude Extents:", longitude_extents)
print("Time Extents:", list(map(dt_to_str, map(_n64_to_datetime, time_extents))))<jupyter_output><empty_output><jupyter_text>Visualize Data Cube Region<jupyter_code>## The code below renders a map that can be used to orient yourself with the region.
from utils.data_cube_utilities.dc_display_map import display_map
display_map(latitude = latitude_extents, longitude = longitude_extents)<jupyter_output><empty_output><jupyter_text>>  Pick a smaller analysis region and display that regionTry to keep your region to less than 0.2-deg x 0.2-deg for rapid processing.Pick the time extents for your mosaic product (keep to 1 year or less).<jupyter_code>## Vietnam - Central Lam Dong Province ##
# longitude_extents = (107.80, 108.00)
# latitude_extents  = (11.70, 11.90)

## Zanzibar - Zanzibar City
latitude_extents  = (-6.25, -6.07)
longitude_extents = (39.15, 39.29)

time_extents = ('2015-01-01', '2015-12-31')
display_map(latitude = latitude_extents, longitude = longitude_extents)<jupyter_output><empty_output><jupyter_text>Load the dataset and the required spectral bands or other parametersAfter loading, you will view the Xarray dataset. Notice the dimensions represent the number of pixels in your latitude and longitude dimension as well as the number of time slices (time) in your time series.<jupyter_code>landsat_dataset = dc.load(latitude = latitude_extents,
                          longitude = longitude_extents,
                          platform = platform,
                          time = time_extents,
                          product = product,
                          measurements = ['red', 'green', 'blue', 'nir', 'swir1', 'swir2', 'pixel_qa']) 
landsat_dataset
#view the dimensions and sample content from the cube<jupyter_output><empty_output><jupyter_text>Display Example Images  > Single band visualization   > For a quick inspection, let's look at two images. The first image will allow the selection of any band (red, blue, green, nir, swir1, swir2) to produce a grey-scale image of any band. The second image will mask clouds with bright red on an RGB image. > Select the desired acquisition (time slice) in the block below. You can select from 1 to , where the max value is the number of time slices noted in the block above. Change the comment statements below to select the bands for the first image.<jupyter_code>acquisition_number = 10
# select an acquisition number from 1 to "time" using the array limits above
%matplotlib inline
# landsat_dataset.red.isel(time = acquisition_number).plot(cmap = "Greys")
landsat_dataset.green.isel(time = acquisition_number).plot(cmap = "Greys")
# landsat_dataset.blue.isel(time = acquisition_number).plot(cmap = "Greys")
#landsat_dataset.nir.isel(time = acquisition_number).plot(cmap = "Greys")
#landsat_dataset.swir1.isel(time = acquisition_number).plot(cmap = "Greys")
#landsat_dataset.swir2.isel(time = acquisition_number).plot(cmap = "Greys")<jupyter_output><empty_output><jupyter_text>> Define Cloud Masking Function   Removes clouds and cloud shadows based on the Landsat pixel QA informationThis is only for reference ... nothing to modify here<jupyter_code>import numpy as np  

def generate_cloud_mask(dataset, include_shadows = False):
    #Create boolean Masks for clear and water pixels
    clear_pixels = dataset.pixel_qa.values == 2  + 64
    water_pixels = dataset.pixel_qa.values == 4  + 64
    shadow_pixels= dataset.pixel_qa.values == 8 + 64
    
    a_clean_mask = np.logical_or(clear_pixels, water_pixels)
    
    if include_shadows:
        a_clean_mask = np.logical_or(a_clean_mask, shadow_pixels)
        
    return np.invert(a_clean_mask)

def remove_clouds(dataset, include_shadows = False):
    #Create boolean Masks for clear and water pixels
    clear_pixels = dataset.pixel_qa.values == 2  + 64
    water_pixels = dataset.pixel_qa.values == 4  + 64
    shadow_pixels= dataset.pixel_qa.values == 8 + 64
    
    a_clean_mask = np.logical_or(clear_pixels, water_pixels)
    
    if include_shadows:
        a_clean_mask = np.logical_or(a_clean_mask, shadow_pixels)
        
    return dataset.where(a_clean_mask)<jupyter_output><empty_output><jupyter_text>>  Mask clouds from your selected acquisition and visualize the scene and mask Now we will look at two RGB images where the second image includes the cloud, cloud shadow and no data mask in RED. Also, the scene is the same as the acquistion selected above.<jupyter_code>cloud_mask = generate_cloud_mask(landsat_dataset)
cloudless = remove_clouds(landsat_dataset)
import matplotlib.pyplot as plt
from utils.data_cube_utilities.dc_rgb import rgb

rgb(landsat_dataset, time_index = acquisition_number)
plt.show()
red = [255,0,0]
rgb(landsat_dataset,time_index = acquisition_number,
    paint_on_mask = [(cloud_mask, red)])<jupyter_output><empty_output><jupyter_text>Cleaning up the clouds and creating a cloud-free mosaicRemember that this process will filter clouds from the entire time series stack >** Most Recent Pixel Mosaic **   >Masks clouds from imagery and uses the most recent cloud-free pixels.<jupyter_code>from utils.data_cube_utilities.dc_mosaic import create_mosaic

def mrf_mosaic(dataset):
    # The mask here is based on pixel_qa products. It comes bundled in with most Landsat Products.
    cloud_free_boolean_mask = np.invert(generate_cloud_mask(dataset))
    return create_mosaic(dataset, clean_mask = cloud_free_boolean_mask)
recent_composite = mrf_mosaic(landsat_dataset)
recent_composite.nir.plot(cmap = "Greys")
rgb(recent_composite)<jupyter_output><empty_output><jupyter_text>----  > **Median Mosaic**  >  Masks clouds from imagery using the median valued cloud-free pixels in the time series<jupyter_code>from utils.data_cube_utilities.dc_mosaic import create_median_mosaic

def median_mosaic(dataset):
    # The mask here is based on pixel_qa products. It comes bundled in with most Landsat Products.
    cloud_free_boolean_mask = np.invert(generate_cloud_mask(dataset))
    return create_median_mosaic(dataset, clean_mask = cloud_free_boolean_mask)
median_composite = median_mosaic(landsat_dataset)
median_composite.nir.plot(cmap = "Greys")
rgb(median_composite)<jupyter_output><empty_output><jupyter_text>Select bands used for clustering<jupyter_code>cluster_bands = ['red', 'green', 'blue', 'swir1']<jupyter_output><empty_output><jupyter_text>Perform K-Means clustering and view the output<jupyter_code>def figure_ratio(ds, fixed_width = 10):
    width = fixed_width
    height = len(ds.latitude) * (fixed_width / len(ds.longitude))
    return (width, height)
from utils.data_cube_utilities.dc_clustering import kmeans_cluster_dataset

# change the number of clusters in the line below, as desired
# this example uses the "median composite" image from above
classification_x =  kmeans_cluster_dataset(median_composite, cluster_bands, n_clusters=8)
# plot the k-mean classification result 
classification_x.plot(figsize = figure_ratio(classification_x))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08575803981623277, 'fraction_numerical': 0.011544351513723643, 'mean_word_length': 5.525749423520369, 'pattern_counts': {'":': 0, '<': 57, '<?xml version=': 0, '>': 72, 'https://': 0, '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': '365537', 'license': 'Apache-2.0', 'n_tokens_mistral': 2787, 'n_tokens_neox': 2645, 'n_words': 735, 'path': 'notebooks/training/ardc_training/Training_TaskA_Mosaics.ipynb', 'repo_name': 'admariner/data_cube_notebooks', 'size': '8489'}