Landsat 8 visualizations
# import geemap
import geemap as geemap
import ee
import numpy as np
import pandas as pd
from egis.utils import (
vis_params_dict, default_vis_params, get_known_roi, extract_single_img,
extract_img_collection_metadata, extract_img_collection_properties,
coords_to_polygon, collection_len
)
ee.Initialize()
Description: understanding Landsat satellite data, frequency of images and ways to visualize without clouds.
Landsat coverage¶
The satellite is circulating around the globe and taking remote sensing images. Thereby it will take consecutive snapshots from the earth as it progresses. Individual areas will be revisited regularly roughly every 2 weeks.
roi_coords = get_known_roi('grand_canyon')
roi_region_large = coords_to_polygon(get_known_roi('us_west_large'))
roi_region_medium = coords_to_polygon(get_known_roi('grand_canyon_rectangle'))
vis_params_l8 = default_vis_params('landsat')
dataset = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR') \
.filterDate('2016-10-01', '2016-12-31') \
.filterBounds(roi_region_large)
The orbit of the satellite is reflected in the image tiles: vertical image neighbours are aligned more than horizontal neighbours, because there is much less time lag between vertically consecutive images.
Map = geemap.Map()
Map.setCenter(roi_coords[0], roi_coords[1], 5)
Map.addLayer(dataset, vis_params_l8)
Map
Due to the time lag, horizontal neighbours can have very different weather conditions.
Map = geemap.Map()
Map.setCenter(roi_coords[0], roi_coords[1], 8)
Map.addLayer(dataset, vis_params_l8)
Map
Image metadata¶
To get a better feeling about the volume of images produced by Landsat we can pull a collection of images and extract metadata from it.
dataset = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR') \
.filterDate('2016-01-01', '2016-12-31') \
.filterBounds(roi_region_medium)
Given those time and spatial filters the collection still has rather many images:
collection_len(dataset)
Any single images comes with a lot of metadata attached:
this_img = extract_single_img(dataset, 0)
this_img_props_dict = geemap.image_props(this_img).getInfo()
img_props = pd.DataFrame.from_dict(this_img_props_dict, orient='index')
img_props
img_props.loc['system:band_names'].squeeze()
We can extract the most important properties for all images:
list_of_props = ['CLOUD_COVER', 'CLOUD_COVER_LAND', 'WRS_PATH', 'WRS_ROW']
all_img_dates = geemap.image_dates(dataset).getInfo()
metadata = extract_img_collection_metadata(dataset, list_of_props)
metadata['date'] = all_img_dates
metadata.head(5)
From this we immediately can see that we get multiple images for multiple tiles of the Worldwide reference system:
metadata.loc[:, ['WRS_PATH', 'WRS_ROW']].drop_duplicates().head(8)
A given tile is revisited roughly every 2 weeks:
single_tile_over_time = metadata.query('WRS_PATH == 39 and WRS_ROW == 37')
single_tile_over_time
single_tile_over_time.shape
Depending on the time of the year, any given tile might look pretty different for different visits of the satellite:
scene = extract_single_img(dataset, 643)
Map = geemap.Map()
Map.centerObject(scene, 8)
Map.addLayer(scene, vis_params_l8, 'default RGB')
Map
In particular with weather phenomena like clouds:
scene = extract_single_img(dataset, 633)
Map = geemap.Map()
Map.centerObject(scene, 8)
Map.addLayer(scene, vis_params_l8, 'default RGB')
Map
Mosaic visualization¶
Clouds are an obvious impediment for any visual inspection. So we can try to improve the visualization of a mosaic by only taking images with rather low ratio of clouds.
Without any special consideration of clouds, this is what we get:
Map = geemap.Map()
Map.setCenter(roi_coords[0], roi_coords[1], 6)
Map.addLayer(dataset, vis_params_l8)
Map
An easy way to improve things is by taking median pixel values, not only the last observation per pixel:
Map = geemap.Map()
Map.setCenter(roi_coords[0], roi_coords[1], 6)
Map.addLayer(dataset.median(), vis_params_l8)
Map
If we sort with regards to cloud coverage, we can easily visualize both extremes for cloud coverage.
cloud_reduced_images = dataset.sort('CLOUD_COVER', opt_ascending=True)
Map = geemap.Map()
Map.setCenter(roi_coords[0], roi_coords[1], 6)
Map.addLayer(cloud_reduced_images, vis_params_l8)
Map
cloud_reduced_images = dataset.sort('CLOUD_COVER', opt_ascending=False)
Map = geemap.Map()
Map.setCenter(roi_coords[0], roi_coords[1], 6)
Map.addLayer(cloud_reduced_images, vis_params_l8)
Map
Although this is already a huge improvement of the visualization, it still comes with the problem that individual images might be taken from totally different points in time. This might lead to neighbouring images having totally different vegetation states (summer vs winter) or even having snow vs no snow.
Masking¶
Masking is a way to create a boolean flag on a pixel level. So for example, for each pixel in an image we could identify whether in this pixel we either see a cloud or a shadow of a cloud. Such a mask can then be used to only show those parts of an image that fulfill a certain level of quality. In the following example we only will show only pixels that are not covered with clouds.
scene = extract_single_img(dataset, 633)
Map = geemap.Map()
Map.centerObject(scene, 8)
Map.addLayer(scene, vis_params_l8, 'default RGB')
Map
def maskL8sr(image):
# From: https://developers.google.com/earth-engine/datasets/catalog/LANDSAT_LC08_C01_T1_SR
# Bits 3 and 5 are cloud shadow and cloud, respectively.
cloudShadowBitMask = (1 << 3)
cloudsBitMask = (1 << 5)
# Get the pixel QA band.
qa = image.select('pixel_qa')
# Both flags should be set to zero, indicating clear conditions.
mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0) \
.And(qa.bitwiseAnd(cloudsBitMask).eq(0))
return image.updateMask(mask)
scene = maskL8sr(extract_single_img(dataset, 633))
Map = geemap.Map()
Map.centerObject(scene, 8)
Map.addLayer(scene, vis_params_l8, 'default RGB')
Map
Now we can also combine masking techniques with median computations: first we filter clouds on a per pixel basis and take the median values of the remaining pixels.
Map = geemap.Map()
Map.setCenter(roi_coords[0], roi_coords[1], 6)
Map.addLayer(dataset.map(maskL8sr).median(), vis_params_l8)
Map
The source of the notebook can be found here