In [1]:
# Enable importing of utilities.
from sys import path
path.append("..")
In [2]:
import datacube
import utils.data_cube_utilities.data_access_api as dc_api
api = dc_api.DataAccessApi()
dc = datacube.Datacube(app = 'ardc_task_a')
api.dc = dc
In [3]:
list_of_products = dc.list_products()
netCDF_products = list_of_products[list_of_products['format'] == 'NetCDF']
netCDF_products
Out[3]:
In [4]:
# Change the data platform and data cube here
platform = "LANDSAT_7"
# platform = "LANDSAT_8"
# product = "ls7_ledaps_ghana"
# product = "ls7_ledaps_kenya"
# product = "ls7_ledaps_senegal"
# product = "ls7_ledaps_sierra_leone"
# product = "ls7_ledaps_tanzania"
product = "ls7_ledaps_vietnam"
In [5]:
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))))
In [6]:
## 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)
Out[6]:
In [7]:
## Vietnam - Central Lam Dong Province ##
longitude_extents = (107.80, 108.00)
latitude_extents = (11.70, 11.90)
time_extents = ('2015-01-01', '2015-12-31')
In [8]:
display_map(latitude = latitude_extents, longitude = longitude_extents)
Out[8]:
In [9]:
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'])
In [10]:
landsat_dataset
#view the dimensions and sample content from the cube
Out[10]:
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.
In [11]:
acquisition_number = 10
# select an acquisition number from 1 to "time" using the array limits above
In [12]:
%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")
Out[12]:
In [13]:
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)
In [14]:
cloud_mask = generate_cloud_mask(landsat_dataset)
cloudless = remove_clouds(landsat_dataset)
In [15]:
import matplotlib.pyplot as plt
from utils.data_cube_utilities.dc_rgb import rgb
rgb(landsat_dataset,at_index = acquisition_number)
plt.show()
In [16]:
red = [255,0,0]
In [17]:
rgb(landsat_dataset,at_index = acquisition_number,paint_on_mask = [(cloud_mask, red)])
Out[17]:
Most Recent Pixel Mosaic
Masks clouds from imagery and uses the most recent cloud-free pixels.
In [18]:
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)
In [19]:
recent_composite = mrf_mosaic(landsat_dataset)
In [20]:
recent_composite.nir.plot(cmap = "Greys")
Out[20]:
In [21]:
rgb(recent_composite)
Out[21]:
Median Mosaic
Masks clouds from imagery using the median valued cloud-free pixels in the time series
In [22]:
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)
In [23]:
median_composite = median_mosaic(landsat_dataset)
In [24]:
median_composite.nir.plot(cmap = "Greys")
Out[24]:
In [25]:
rgb(median_composite)
Out[25]:
In [26]:
cluster_bands = ['red', 'green', 'blue', 'swir1']
In [27]:
def figure_ratio(ds, fixed_width = 10):
width = fixed_width
height = len(ds.latitude) * (fixed_width / len(ds.longitude))
return (width, height)
In [28]:
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)
In [29]:
# plot the k-mean classification result
classification_x.plot(figsize = figure_ratio(classification_x))
Out[29]: