In [1]:
import matplotlib.pyplot as plt
%matplotlib inline
# supress warnings
import warnings
warnings.filterwarnings('ignore')
In [2]:
def NDVI(dataset):
return (dataset.nir - dataset.red)/(dataset.nir + dataset.red)
In [3]:
def NDWI(dataset):
return (dataset.green - dataset.nir)/(dataset.green + dataset.nir)
In [4]:
def NDBI(dataset):
return (dataset.swir2 - dataset.nir)/(dataset.swir2 + dataset.nir)
In [5]:
from utils.data_cube_utilities.dc_mosaic import create_mosaic, ls8_unpack_qa
import numpy as np
def mosaic(dataset):
# The mask here is based on pixel_qa. It comes bundled in with most Landsat Products.
clear_xarray = ls8_unpack_qa(dataset.pixel_qa, "clear") # Boolean Xarray indicating landcover
water_xarray = ls8_unpack_qa(dataset.pixel_qa, "water") # Boolean Xarray indicating watercover
cloud_free_boolean_mask = np.logical_or(clear_xarray, water_xarray)
return create_mosaic(dataset, clean_mask = cloud_free_boolean_mask)
Median Mosaic
A cloud free representation fo satellite imagery. Works by masking out clouds from imagery, and using the median valued cloud-free pixels in the time series
In [6]:
from utils.data_cube_utilities.dc_mosaic import create_median_mosaic, ls8_unpack_qa
def median_mosaic(dataset):
# The mask here is based on pixel_qa. It comes bundled in with most Landsat Products.
clear_xarray = ls8_unpack_qa(dataset.pixel_qa, "clear") # Boolean Xarray indicating landcover
water_xarray = ls8_unpack_qa(dataset.pixel_qa, "water") # Boolean Xarray indicating watercover
cloud_free_boolean_mask = np.logical_or(clear_xarray, water_xarray)
return create_median_mosaic(dataset, clean_mask = cloud_free_boolean_mask)
Data cube object
A datacube object is your interface with data stored on your data cube system.
In [7]:
import datacube
dc = datacube.Datacube(app = '3B_urban')
Loading a Dataset
Requires latitude-longitude bounds of an area, a time-range, list of desired measurements, platform and product names.
In [8]:
dc.list_products()
Out[8]:
In [9]:
## Loading in a region
start_date = "2019-01-01"
end_date = "2019-12-31"
# Kumasi, Ghana
lat = (6.597724,6.781856)
lon = (-1.727843,-1.509147)
## Considering the year of 2015
date_range = (start_date,end_date)
## Landsat 8 Data
platform = 'LANDSAT_8'
product = 'ls8_usgs_sr_scene'
desired_bands = ['red','green','nir','swir2', 'pixel_qa'] # needed by ndvi, ndwi, ndbi and cloud masking
desired_bands = desired_bands + ['blue'] # blue is needed for a true color visualization purposes
landsat_dataset = dc.load(product = product,\
platform = platform,\
lat = lat,\
lon = lon,\
time = date_range,\
measurements = desired_bands)
In [10]:
landsat_mosaic = median_mosaic(landsat_dataset)
Saving your data
A .tiff or png is a great way to represent true color mosaics. The image below is a saved .png representation of of a landsat mosaic.
In [11]:
landsat_mosaic
Out[11]:
In [12]:
import os
import pathlib
from utils.data_cube_utilities.dc_utilities import write_png_from_xr
png_dir = 'output/pngs'
pathlib.Path(png_dir).mkdir(parents=True, exist_ok=True)
write_png_from_xr(png_dir + '/cloud_free_mosaic.png',
landsat_mosaic, ["red", "green", "blue"],
scale = [(0,2000),(0,2000),(0,2000)])
In [13]:
ndbi = NDBI(landsat_mosaic) # Urbanization
ndvi = NDVI(landsat_mosaic) # Dense Vegetation
ndwi = NDWI(landsat_mosaic) # High Concentrations of Water
Plot Values
xarray data-arrays have built in plotting functions you can use to validate trends or differences in your data.
In [14]:
ndvi.plot(cmap = "Greens")
plt.show()
In [15]:
ndwi.plot(cmap = "Blues")
plt.show()
In [16]:
(ndbi + 0.2).plot(cmap = "Reds")
plt.show()
Convert To a Dataset
It's good practice to accurately name your datasets and data-arrays. If you'd like to merge data-arrays into a larger datasets, you should convert data-arrays to datasets
In [17]:
ds_ndvi = ndvi.to_dataset(name = "NDVI")
ds_ndwi = ndwi.to_dataset(name= "NDWI")
ds_ndbi = ndbi.to_dataset(name = "NDBI")
Merge into one large Dataset
If your data-arrays share the same set of coordinates, or if you feel that you'll be using these values together in the future, you should consider merging them into a dataset
In [18]:
urbanization_dataset = ds_ndvi.merge(ds_ndwi).merge(ds_ndbi)
Checking your Merge
The string readout of your new dataset should give you a good idea about how the.mergewent.
In [19]:
urbanization_dataset
Out[19]:
Building a False Color Composite
If you have three lowly correlated measurements, place each measurement on its own Red, Green, Blue channel and visualize it.
In [20]:
write_png_from_xr(png_dir + '/false_color.png', urbanization_dataset,
["NDBI", "NDVI", "NDWI"], scale = [(-1,1),(0,1),(0,1)])
Analyze The False Color Image
Values that adhere strongly to individual classes adhere to their own color channel. In this example, NDVI adheres to green, NDWI adheres to blue, and NDBI seems to adhere to red
Validate urbanization using other imagery
Double check results using high-resolution imagery. Compare to the false color mosaic
In [21]:
from utils.data_cube_utilities.dc_display_map import display_map
display_map(latitude = lat ,longitude = lon)
Out[21]: