This notebook can be used to create custom Landsat cloud-filtered mosaics for any time period and location. The mosaics can be output as GeoTIFF products for analysis in external GIS tools. The following mosaics are possible:
Median = midpoint of spectral data Geomedian = Australian median product with improved spectral consistency Most-Recent = most-recent clear pixel Max-NDVI = maximum vegetation response
Users should review the "Cloud_Statistics" notebook for more information about the cloud statistics for any given temporal and spatial combination. An understanding of the underlying data is important for creating a valid mosaic for further analyses.
In [1]:
# Enable importing of utilities.
import sys
sys.path.append('..')
# Supress Warnings.
import warnings
warnings.filterwarnings('ignore')
import numpy as np
import xarray as xr
import pandas as pd
import matplotlib.pyplot as plt
# Load Data Cube Configuration
import datacube
import utils.data_cube_utilities.data_access_api as dc_api
api = dc_api.DataAccessApi()
dc = api.dc
List available products for each platform
In [2]:
# Get available products
products_info = dc.list_products()
# List LANDSAT 7 products
print("LANDSAT 7 Products:")
products_info[["platform", "name"]][products_info.platform == "LANDSAT_7"]
Out[2]:
In [3]:
# List LANDSAT 8 products
print("LANDSAT 8 Products:")
products_info[["platform", "name"]][products_info.platform == "LANDSAT_8"]
Out[3]:
Choose product
<p style="color:red";>CHANGE INPUTS BELOW
In [5]:
# Select a Product and Platform
# Examples: ghana, kenya, tanzania, sierra_leone, senegal
# product = "ls7_usgs_sr_scene"
# platform = "LANDSAT_7"
product = "ls8_usgs_sr_scene"
platform = "LANDSAT_8"
In [6]:
from utils.data_cube_utilities.dc_load import get_product_extents
from utils.data_cube_utilities.dc_time import dt_to_str
full_lat, full_lon, min_max_dates = get_product_extents(api, platform, product)
# Print the extents of the combined data.
print("Latitude Extents:", full_lat)
print("Longitude Extents:", full_lon)
print("Time Extents:", list(map(dt_to_str, min_max_dates)))
Visualize the available area
In [7]:
from utils.data_cube_utilities.dc_display_map import display_map
display_map(full_lat, full_lon)
Out[7]:
<p style="color:red";>CHANGE INPUTS BELOW
In [14]:
# Select an analysis region (Lat-Lon) within the extents listed above.
# HINT: Keep your region small (<0.5 deg square) to avoid memory overload issues
# Select a time period (Min-Max) within the extents listed above (Year-Month-Day)
# This region and time period will be used for the cloud assessment
# Mombasa, Kenya
latitude = (-4.1, -3.9)
longitude = (39.5, 39.7)
# Freetown, Sierra Leone
# latitude = (8.3267, 8.5123)
# longitude = (-13.3109, -13.1197 )
# Tano-Offin Forest - Ghana
# latitude = (6.5814, 6.8978 )
# longitude = (-2.2955, -1.9395)
# Vietnam
# latitude = (10.9358, 11.0358)
# longitude = (107.1899, 107.2899)
# Time Period
time_extents = ('2018-01-01', '2018-12-31')
Visualize the selected area
In [15]:
display_map(latitude,longitude)
Out[15]:
In [16]:
landsat_dataset = dc.load(latitude = latitude,
longitude = longitude,
platform = platform,
time = time_extents,
product = product,
measurements = ['red', 'green', 'blue', 'nir', 'swir1', 'swir2', 'pixel_qa'])
In [17]:
# Displays the first few values of each data array to check the content
# Latitude and Longitude numbers = number of pixels in each dimension
# Time = number of time slices in the dataset
landsat_dataset
Out[17]:
In [18]:
from utils.data_cube_utilities.clean_mask import landsat_qa_clean_mask
cloud_mask = landsat_qa_clean_mask(landsat_dataset, platform=platform)
cleaned_dataset = landsat_dataset.where(cloud_mask)
Median Mosaic
Masks clouds from imagery using the median-valued cloud-free pixels in the time series. More specifically, each band (e.g. red) of each pixel is assigned its median across time. So this mosaic method generates values that are not in the dataset.
In [19]:
from utils.data_cube_utilities.dc_mosaic import create_median_mosaic
from utils.data_cube_utilities.dc_rgb import rgb
In [20]:
# RGB image options
# Standard RGB = 321 = Red, Green, Blue
# False Color = 543 = SWIR1, NIR, Red
# False Color (Landsat Mosaic) = 742 = SWIR2, NIR, Green
median_composite = create_median_mosaic(cleaned_dataset, cloud_mask)
rgb(median_composite, bands=['swir2', 'nir', 'green'], min_possible=0, max_possible=4000, width=10)
plt.show()
Geomedian Mosaic
Masks clouds from imagery using the geomedian-valued cloud-free pixels in the time series, which maintains the spectral band relationships. That is, this is a median through time for all bands considered collectively rather than separately, as is the case in a median mosaic. This algorithm was developed by Geoscience Australia and produces excellent cloud-filtered mosaic products for further analysis.
For more information, see the following paper: High-Dimensional Pixel Composites from Earth Observation Time Series, by, Dale Roberts, Norman Mueller, and Alexis McIntyre. IEEE Transactions on Geoscience and Remote Sensing, Vol. 55. No. 11, November 2017.
In [21]:
from utils.data_cube_utilities.dc_mosaic import create_hdmedians_multiple_band_mosaic
In [22]:
# RGB image options
# Standard RGB = 321 = Red, Green, Blue
# False Color = 543 = SWIR1, NIR, Red
# False Color (Landsat Mosaic) = 742 = SWIR2, NIR, Green
geomedian_composite = create_median_mosaic(cleaned_dataset, cloud_mask)
rgb(geomedian_composite, bands=['swir2', 'nir', 'green'], min_possible=0, max_possible=4000, width=10)
plt.show()
Most Recent and Least Recent Mosaic
Masks clouds from imagery using the most or least recent cloud-free pixels in the time series.
In [23]:
from utils.data_cube_utilities.dc_mosaic import create_mosaic
Most Recent Mosaic
In [24]:
# RGB image options
# Standard RGB = 321 = Red, Green, Blue
# False Color = 543 = SWIR1, NIR, Red
# False Color (Landsat Mosaic) = 742 = SWIR2, NIR, Green
most_recent_composite = create_mosaic(cleaned_dataset, cloud_mask.values)
rgb(most_recent_composite, bands=['swir2', 'nir', 'green'], min_possible=0, max_possible=4000, width=10)
plt.show()
Least Recent Mosaic
In [25]:
least_recent_composite = create_mosaic(cleaned_dataset, cloud_mask.values, reverse_time=True)
rgb(least_recent_composite, bands=['swir2', 'nir', 'green'], min_possible=0, max_possible=4000, width=10)
plt.show()
Max NDVI Mosaic
Masks clouds from imagery using the Max NDVI across time for cloud-free pixels in the time series.
In [26]:
from utils.data_cube_utilities.dc_mosaic import create_max_ndvi_mosaic
In [27]:
# RGB image options
# Standard RGB = 321 = Red, Green, Blue
# False Color = 543 = SWIR1, NIR, Red
# False Color (Landsat Mosaic) = 742 = SWIR2, NIR, Green
max_ndvi_composite = create_max_ndvi_mosaic(cleaned_dataset, cloud_mask.values)
rgb(max_ndvi_composite, bands=['swir2', 'nir', 'green'], min_possible=0, max_possible=4000, width=10)
plt.show()
In [28]:
import os
from utils.data_cube_utilities.import_export import export_slice_to_geotiff
# Remove the comment tags (#) to export a GeoTIFF output product
# Change the name of the output file, or it will be overwritten for each run
output_dir = 'output/geotiffs'
if not os.path.exists(output_dir):
os.makedirs(output_dir)
export_slice_to_geotiff(median_composite, output_dir + '/DEMO_median_composite.tif')
# export_slice_to_geotiff(geomedian_composite, output_dir + '/DEMO_geomedian_composite.tif')
# export_slice_to_geotiff(most_recent_composite, output_dir + '/DEMO_most_recent_composite.tif')
# export_slice_to_geotiff(max_ndvi_composite, output_dir + '/DEMO_max_ndvi_composite.tif')
In [29]:
!ls -lah output/geotiffs/