Spectral Products: Vegetation, Water, Urbanization¶
In [1]:
# Supress Warning
import warnings
warnings.filterwarnings('ignore')
# Load Data Cube Configuration
import datacube
dc = datacube.Datacube(app = 'my_app')
import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np
import xarray as xr
# Import Data Cube API
import utils.data_cube_utilities.data_access_api as dc_api
api = dc_api.DataAccessApi()
In [2]:
## LS8 Caqueta
# Latitude: (0.000134747292617865, 1.077843593651382)
# Longitude: (-74.91935994831539, -73.30266193148462)
# '2013-04-13', '2018-03-26'
# Resolution: (-0.000269494585236, 0.000269494585236)
## LS8 Vietnam
# Latitude: (10.513927001104687, 12.611133863411238)
# Longitude: (106.79005909290998, 108.91906631627438)
# '2014-01-14', '2016-12-21'
# Resolution: (-0.000269494585236, 0.000269494585236)
## LS7 Caqueta
# Latitude: (0.000134747292617865, 1.077843593651382)
# Longitude: (-74.91935994831539, -73.30266193148462)
# '1999-08-21', '2018-03-25'
# Resolution: (-0.000269494585236, 0.000269494585236)
## LS7 Lake Baringo
# Latitude: (0.4997747685, 0.7495947795)
# Longitude: (35.9742163305, 36.473586859499996)
# '2005-01-08', '2016-12-24'
# Resolution: (-0.000269493, 0.000269493)
In [3]:
# CHANGE HERE >>>>>>>>>>>>>>>>>
# Select a Product and Platform
# product = ls7_collection1_AMA_ingest
# platform = "LANDSAT_7"
product = 'ls8_collection1_AMA_ingest'
platform = "LANDSAT_8"
output_crs = 'EPSG:4326'
resolution = (-0.000269494585236, 0.000269494585236)
In [4]:
# CHANGE HERE >>>>>>>>>>>>>>>>>>
# 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
#Sub-region selection
# LS7 Lake Baringo
#latitude = (0.49964002, 0.76)
#longitude = (36.0, 36.16)
#time_extents = ('2015-01-01', '2018-01-01')
#LS8 Vietnam
latitude = (11.3124, 10.9124)
longitude = (106.9170, 107.3170)
time_extents = ('2014-01-01', '2016-01-01')
In [5]:
# The code below renders a map that can be used to view the region.
from utils.data_cube_utilities.dc_display_map import display_map
display_map(latitude,longitude)
Out[5]:
Load the dataset and the required spectral bands or other parameters¶
In [6]:
landsat_dataset = dc.load(latitude = latitude,
longitude = longitude,
platform = platform,
time = time_extents,
product = product,
output_crs=output_crs,
resolution=resolution,
measurements = ['red', 'green', 'blue', 'nir', 'swir1', 'swir2', 'pixel_qa'])
In [7]:
# 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[7]:
Mask out clouds and cloud shadows + water (if desired) and create a median mosaic¶
In [8]:
from utils.data_cube_utilities.clean_mask import landsat_qa_clean_mask
cloud_mask = landsat_qa_clean_mask(landsat_dataset, platform=platform)
land_mask = landsat_qa_clean_mask(landsat_dataset, platform=platform, cover_types=['clear'])
# Land and Water Dataset = Land and Water pixels with NO Clouds and NO Cloud Shadows
land_and_water_dataset = landsat_dataset.where(cloud_mask)
# Land Dataset = Land ONLY pixels with NO Clouds, NO Cloud Shadows and NO Water pixels
land_dataset = landsat_dataset.where(land_mask)
from utils.data_cube_utilities.dc_mosaic import create_median_mosaic, create_max_ndvi_mosaic, create_hdmedians_multiple_band_mosaic
In [9]:
# CHANGE HERE >>>>>>>>>>>>>>>>>>>>
# Select a compositing method to create your cloud-filtered mosaic
# Remove the comments from the pair of lines under one of the mosaic types
# Options are: Median, Geomedian, or Max_NDVI
# This is the MEDIAN mosaic
land_and_water_composite = create_median_mosaic(land_and_water_dataset, cloud_mask)
land_composite = create_median_mosaic(land_dataset, land_mask)
# This is the GEOMEDIAN mosaic
# land_and_water_composite = create_hdmedians_multiple_band_mosaic(land_and_water_dataset, cloud_mask)
# land_composite = create_hdmedians_multiple_band_mosaic(land_dataset, land_mask)
# This is the MAX_NDVI mosaic
# land_and_water_composite = create_max_ndvi_mosaic(land_and_water_dataset, cloud_mask)
# land_composite = create_max_ndvi_mosaic(land_dataset, land_mask)
In [10]:
# Show the land and water composite
from utils.data_cube_utilities.dc_rgb import rgb
rgb(land_and_water_composite)
Out[10]:
In [11]:
# Show the land-only composite (water will be removed ... black pixels)
from utils.data_cube_utilities.dc_rgb import rgb
rgb(land_composite)
Out[11]:
Land Fractional Cover¶
Fractional Cover (FC) is used for landcover type estimation (vegetation, non-green vegetation, bare soil) of each pixel.
We use a model from CSIRO (Juan Gerschmann) and apply it to a median mosaic.
Bare Soil (bs), Photosynthetic Vegetation (pv), Non Photosynthetic Vegetation (npv)
Plot a False Color RGB result where RGB = bs/pv/npv
In [12]:
from utils.data_cube_utilities.dc_fractional_coverage_classifier import frac_coverage_classify
frac_classes = frac_coverage_classify(land_composite, clean_mask = np.ones(land_composite.pixel_qa.shape).astype(np.bool))
In [13]:
# Plot of Fractional Cover
# RED = Bare Soil or Urban Areas
# BLUE = Non-Green Vegetation
# GREEN = Green Vegetation
# BLACK = Water
rgb(frac_classes, bands = ['bs', 'pv', 'npv'], use_data_min=True, use_data_max=True)
Out[13]: