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_b')
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]:
Pick a smaller analysis region and display that region
Try to keep your region to less than 0.2-deg x 0.2-deg for rapid processing. You can click on the map above to find the Lat-Lon coordinates of any location. You will want to identify a region with an inland water body. Pick a time window of a few months so we can pick out some clear pixels and plot the water.
In [8]:
## Vietnam - Central Lam Dong Province ##
longitude_extents = (107.0, 107.2)
latitude_extents = (11.7, 12.0)
time_extents = ('2010-01-01', '2010-12-31')
In [9]:
display_map(latitude = latitude_extents, longitude = longitude_extents)
Out[9]:
In [10]:
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 [11]:
landsat_dataset
#view the dimensions and sample content from the cube
Out[11]:
Single band visualization
For a quick inspection, let's look at one image. The code will allow the selection of any band (red, blue, green, nir, swir1, swir2) to produce a grey-scale 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 [12]:
acquisition_number = 2
# select an acquisition number from 1 to "time" using the array limits above
In [13]:
%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[13]:
In [14]:
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 [15]:
cloud_mask = generate_cloud_mask(landsat_dataset)
cloudless = remove_clouds(landsat_dataset) #landsat_dataset.where(image_is_clean)
Set up plotting function (to be used later) Nothing to modify here
In [16]:
from utils.data_cube_utilities.dc_rgb import rgb
Most Recent Pixel Mosaic
Masks clouds from imagery and uses the most recent cloud-free pixels.
In [17]:
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 [18]:
recent_composite = mrf_mosaic(landsat_dataset)
In [19]:
recent_composite.nir.plot(cmap = "Greys")
Out[19]:
In [20]:
rgb(recent_composite, width = 20)
Out[20]:
In [21]:
from utils.data_cube_utilities.dc_water_classifier import wofs_classify
In [22]:
water_classification = wofs_classify(recent_composite, clean_mask = np.ones(recent_composite.pixel_qa.shape).astype(np.bool), mosaic = True)
In [23]:
water_classification.wofs.plot(cmap='Blues')
Out[23]:
In [24]:
def NDWI(dataset):
return (dataset.green - dataset.nir)/(dataset.green + dataset.nir)
In [25]:
ndwi = NDWI(recent_composite) # High Concentrations of Water - Blues
In [26]:
(ndwi).plot(cmap = "Blues")
Out[26]:
In [27]:
from utils.data_cube_utilities.dc_water_quality import tsm
In [28]:
mask_that_only_includes_water_pixels = water_classification.wofs == 1
tsm_dataset = tsm(recent_composite, clean_mask = mask_that_only_includes_water_pixels )
In [29]:
tsm_dataset.tsm.plot(cmap = "jet")
Out[29]:
In [30]:
ts_water_classification = wofs_classify(landsat_dataset,clean_mask = np.invert(cloud_mask))
In [31]:
# Apply nan to no_data values
ts_water_classification = ts_water_classification.where(ts_water_classification != -9999)
##Time series aggregation that ignores nan values.
water_classification_percentages = (ts_water_classification.mean(dim = ['time']) * 100).wofs.rename('water_classification_percentages')
In [32]:
## import color-scheme and set nans to black
from matplotlib.cm import jet_r as jet_r
jet_r.set_bad('black',1)
## apply nan to percentage values that aren't greater than 0, then plot
water_classification_percentages\
.where(water_classification_percentages > 0)\
.plot(cmap = jet_r)
Out[32]:
In [33]:
pixel_lat = 11.84
pixel_lon = 107.09
In [34]:
pixel = ts_water_classification.sel( latitude = pixel_lat,
longitude = pixel_lon,
method = 'nearest') # nearest neighbor selection
In [35]:
import matplotlib.pyplot as plt
plt.figure(figsize = (20,5))
plt.scatter(pixel.time.values, pixel.wofs.values)
Out[35]: