## dea_plotting.py
'''
Plotting and animating Digital Earth Australia products and data.
License: The code in this notebook is licensed under the Apache License,
Version 2.0 (https://www.apache.org/licenses/LICENSE-2.0). Digital Earth
Australia data is licensed under the Creative Commons by Attribution 4.0
license (https://creativecommons.org/licenses/by/4.0/).
Contact: If you need assistance, please post a question on the Open Data
Cube Slack channel (http://slack.opendatacube.org/) or on the GIS Stack
Exchange (https://gis.stackexchange.com/questions/ask?tags=open-data-cube)
using the `open-data-cube` tag (you can view previously asked questions
here: https://gis.stackexchange.com/questions/tagged/open-data-cube).
If you would like to report an issue with this script, file one on
GitHub: https://github.com/GeoscienceAustralia/dea-notebooks/issues/new
Last modified: April 2023
'''
# Import required packages
import math
import folium
import numpy as np
import pandas as pd
import geopandas as gpd
import matplotlib.patheffects as PathEffects
import matplotlib.pyplot as plt
import xarray as xr
from matplotlib import colors as mcolours
from matplotlib.animation import FuncAnimation
from pathlib import Path
from pyproj import Transformer
from shapely.geometry import box
from skimage.exposure import rescale_intensity
from tqdm.auto import tqdm
import odc.geo.xr
from odc.ui import image_aspect
from dea_tools.spatial import add_geobox
[docs]
def rgb(ds,
bands=['nbart_red', 'nbart_green', 'nbart_blue'],
index=None,
index_dim='time',
robust=True,
percentile_stretch=None,
col_wrap=4,
size=6,
aspect=None,
titles=None,
savefig_path=None,
savefig_kwargs={},
**kwargs):
"""
Takes an xarray dataset and plots RGB images using three imagery
bands (e.g ['nbart_red', 'nbart_green', 'nbart_blue']). The `index`
parameter allows easily selecting individual or multiple images for
RGB plotting. Images can be saved to file by specifying an output
path using `savefig_path`.
This function was designed to work as an easier-to-use wrapper
around xarray's `.plot.imshow()` functionality.
Last modified: September 2020
Parameters
----------
ds : xarray Dataset
A two-dimensional or multi-dimensional array to plot as an RGB
image. If the array has more than two dimensions (e.g. multiple
observations along a 'time' dimension), either use `index` to
select one (`index=0`) or multiple observations
(`index=[0, 1]`), or create a custom faceted plot using e.g.
`col="time"`.
bands : list of strings, optional
A list of three strings giving the band names to plot. Defaults
to '['nbart_red', 'nbart_green', 'nbart_blue']'.
index : integer or list of integers, optional
`index` can be used to select one (`index=0`) or multiple
observations (`index=[0, 1]`) from the input dataset for
plotting. If multiple images are requested these will be plotted
as a faceted plot.
index_dim : string, optional
The dimension along which observations should be plotted if
multiple observations are requested using `index`. Defaults to
`time`.
robust : bool, optional
Produces an enhanced image where the colormap range is computed
with 2nd and 98th percentiles instead of the extreme values.
Defaults to True.
percentile_stretch : tuple of floats
An tuple of two floats (between 0.00 and 1.00) that can be used
to clip the colormap range to manually specified percentiles to
get more control over the brightness and contrast of the image.
The default is None; '(0.02, 0.98)' is equivelent to
`robust=True`. If this parameter is used, `robust` will have no
effect.
col_wrap : integer, optional
The number of columns allowed in faceted plots. Defaults to 4.
size : integer, optional
The height (in inches) of each plot. Defaults to 6.
aspect : integer, optional
Aspect ratio of each facet in the plot, so that aspect * size
gives width of each facet in inches. Defaults to None, which
will calculate the aspect based on the x and y dimensions of
the input data.
titles : string or list of strings, optional
Replace the xarray 'time' dimension on plot titles with a string
or list of string titles, when a list of index values are
provided, of your choice. Defaults to None.
savefig_path : string, optional
Path to export image file for the RGB plot. Defaults to None,
which does not export an image file.
savefig_kwargs : dict, optional
A dict of keyword arguments to pass to
`matplotlib.pyplot.savefig` when exporting an image file. For
all available options, see:
https://matplotlib.org/api/_as_gen/matplotlib.pyplot.savefig.html
**kwargs : optional
Additional keyword arguments to pass to `xarray.plot.imshow()`.
For example, the function can be used to plot into an existing
matplotlib axes object by passing an `ax` keyword argument.
For more options, see:
http://xarray.pydata.org/en/stable/generated/xarray.plot.imshow.html
Returns
-------
An RGB plot of one or multiple observations, and optionally an image
file written to file.
"""
# Get names of x and y dims
y_dim, x_dim = ds.odc.spatial_dims
# If ax is supplied via kwargs, ignore aspect and size
if 'ax' in kwargs:
# Create empty aspect size kwarg that will be passed to imshow
aspect_size_kwarg = {}
else:
# Compute image aspect
if not aspect:
aspect = ds.odc.geobox.aspect
# Populate aspect size kwarg with aspect and size data
aspect_size_kwarg = {'aspect': aspect, 'size': size}
# If no value is supplied for `index` (the default), plot using default
# values and arguments passed via `**kwargs`
if index is None:
# Select bands and convert to DataArray
da = ds[bands].to_array().compute()
# If percentile_stretch == True, clip plotting to percentile vmin, vmax
if percentile_stretch:
vmin, vmax = da.quantile(percentile_stretch).values
kwargs.update({'vmin': vmin, 'vmax': vmax})
# If there are more than three dimensions and the index dimension == 1,
# squeeze this dimension out to remove it
if ((len(ds.dims) > 2) and ('col' not in kwargs) and
(len(da[index_dim]) == 1)):
da = da.squeeze(dim=index_dim)
# If there are more than three dimensions and the index dimension
# is longer than 1, raise exception to tell user to use 'col'/`index`
elif ((len(ds.dims) > 2) and ('col' not in kwargs) and
(len(da[index_dim]) > 1)):
raise Exception(
f'The input dataset `ds` has more than two dimensions: '
f'{list(ds.dims.keys())}. Please select a single observation '
'using e.g. `index=0`, or enable faceted plotting by adding '
'the arguments e.g. `col="time", col_wrap=4` to the function '
'call')
img = da.plot.imshow(x=x_dim,
y=y_dim,
robust=robust,
col_wrap=col_wrap,
**aspect_size_kwarg,
**kwargs)
if titles is not None:
for ax, title in zip(img.axs.flat, titles):
ax.set_title(title)
# If values provided for `index`, extract corresponding observations and
# plot as either single image or facet plot
else:
# If a float is supplied instead of an integer index, raise exception
if isinstance(index, float):
raise Exception(
f'Please supply `index` as either an integer or a list of '
'integers')
# If col argument is supplied as well as `index`, raise exception
if 'col' in kwargs:
raise Exception(
f'Cannot supply both `index` and `col`; please remove one and '
'try again')
# Convert index to generic type list so that number of indices supplied
# can be computed
index = index if isinstance(index, list) else [index]
# Select bands and observations and convert to DataArray
da = ds[bands].isel(**{index_dim: index}).to_array().compute()
# If percentile_stretch == True, clip plotting to percentile vmin, vmax
if percentile_stretch:
vmin, vmax = da.quantile(percentile_stretch).values
kwargs.update({'vmin': vmin, 'vmax': vmax})
# If multiple index values are supplied, plot as a faceted plot
if len(index) > 1:
img = da.plot.imshow(x=x_dim,
y=y_dim,
robust=robust,
col=index_dim,
col_wrap=col_wrap,
**aspect_size_kwarg,
**kwargs)
if titles is not None:
for ax, title in zip(img.axs.flat, titles):
ax.set_title(title)
# If only one index is supplied, squeeze out index_dim and plot as a
# single panel
else:
img = da.squeeze(dim=index_dim).plot.imshow(robust=robust,
**aspect_size_kwarg,
**kwargs)
if titles is not None:
for ax, title in zip(img.axs.flat, titles):
ax.set_title(title)
# If an export path is provided, save image to file. Individual and
# faceted plots have a different API (figure vs fig) so we get around this
# using a try statement:
if savefig_path:
print(f'Exporting image to {savefig_path}')
try:
img.fig.savefig(savefig_path, **savefig_kwargs)
except:
img.figure.savefig(savefig_path, **savefig_kwargs)
[docs]
def display_map(x, y, crs='EPSG:4326', margin=-0.5, zoom_bias=0):
"""
Given a set of x and y coordinates, this function generates an
interactive map with a bounded rectangle overlayed on Google Maps
imagery.
Last modified: February 2023
Modified from function written by Otto Wagner available here:
https://github.com/ceos-seo/data_cube_utilities/tree/master/data_cube_utilities
Parameters
----------
x : (float, float)
A tuple of x coordinates in (min, max) format.
y : (float, float)
A tuple of y coordinates in (min, max) format.
crs : string, optional
A string giving the EPSG CRS code of the supplied coordinates.
The default is 'EPSG:4326'.
margin : float
A numeric value giving the number of degrees lat-long to pad
the edges of the rectangular overlay polygon. A larger value
results more space between the edge of the plot and the sides
of the polygon. Defaults to -0.5.
zoom_bias : float or int
A numeric value allowing you to increase or decrease the zoom
level by one step. Defaults to 0; set to greater than 0 to zoom
in, and less than 0 to zoom out.
Returns
-------
folium.Map : A map centered on the supplied coordinate bounds. A
rectangle is drawn on this map detailing the perimeter of the x, y
bounds. A zoom level is calculated such that the resulting
viewport is the closest it can possibly get to the centered
bounding rectangle without clipping it.
"""
# Convert each corner coordinates to lat-lon
all_x = (x[0], x[1], x[0], x[1])
all_y = (y[0], y[0], y[1], y[1])
transformer = Transformer.from_crs(crs, "EPSG:4326")
all_longitude, all_latitude = transformer.transform(all_x, all_y)
# Calculate zoom level based on coordinates
lat_zoom_level = _degree_to_zoom_level(
min(all_latitude), max(all_latitude), margin=margin) + zoom_bias
lon_zoom_level = _degree_to_zoom_level(
min(all_longitude), max(all_longitude), margin=margin) + zoom_bias
zoom_level = min(lat_zoom_level, lon_zoom_level)
# Identify centre point for plotting
center = [np.mean(all_latitude), np.mean(all_longitude)]
# Create map
interactive_map = folium.Map(
location=center,
zoom_start=zoom_level,
tiles="http://mt1.google.com/vt/lyrs=y&z={z}&x={x}&y={y}",
attr="Google")
# Create bounding box coordinates to overlay on map
line_segments = [(all_latitude[0], all_longitude[0]),
(all_latitude[1], all_longitude[1]),
(all_latitude[3], all_longitude[3]),
(all_latitude[2], all_longitude[2]),
(all_latitude[0], all_longitude[0])]
# Add bounding box as an overlay
interactive_map.add_child(
folium.features.PolyLine(locations=line_segments,
color='red',
opacity=0.8))
# Add clickable lat-lon popup box
interactive_map.add_child(folium.features.LatLngPopup())
return interactive_map
[docs]
def xr_animation(ds,
bands=None,
output_path='animation.mp4',
width_pixels=500,
interval=100,
percentile_stretch=(0.02, 0.98),
image_proc_funcs=None,
show_gdf=None,
show_date='%d %b %Y',
show_text=None,
show_colorbar=True,
gdf_kwargs={},
annotation_kwargs={},
imshow_kwargs={},
colorbar_kwargs={},
limit=None):
"""
Takes an `xarray` timeseries and animates the data as either a
three-band (e.g. true or false colour) or single-band animation,
allowing changes in the landscape to be compared across time.
Animations can be customised to include text and date annotations
or use specific combinations of input bands. Vector data can be
overlaid and animated on top of imagery, and custom image
processing functions can be applied to each frame.
Supports .mp4 (ideal for Twitter/social media) and .gif (ideal
for all purposes, but can have large file sizes) format files.
Last modified: April 2023
Parameters
----------
ds : xarray.Dataset
An xarray dataset with multiple time steps (i.e. multiple
observations along the `time` dimension).
bands : list of strings
An list of either one or three band names to be plotted,
all of which must exist in `ds`.
output_path : str, optional
A string giving the output location and filename of the
resulting animation. File extensions of '.mp4' and '.gif' are
accepted. Defaults to 'animation.mp4'.
width_pixels : int, optional
An integer defining the output width in pixels for the
resulting animation. The height of the animation is set
automatically based on the dimensions/ratio of the input
xarray dataset. Defaults to 500 pixels wide.
interval : int, optional
An integer defining the milliseconds between each animation
frame used to control the speed of the output animation. Higher
values result in a slower animation. Defaults to 100
milliseconds between each frame.
percentile_stretch : tuple of floats, optional
An optional tuple of two floats that can be used to clip one or
three-band arrays by percentiles to produce a more vibrant,
visually attractive image that is not affected by outliers/
extreme values. The default is `(0.02, 0.98)` which is
equivalent to xarray's `robust=True`. This parameter is ignored
completely if `vmin` and `vmax` are provided as kwargs to
`imshow_kwargs`.
image_proc_funcs : list of funcs, optional
An optional list containing functions that will be applied to
each animation frame (timestep) prior to animating. This can
include image processing functions such as increasing contrast,
unsharp masking, saturation etc. The function should take AND
return a `numpy.ndarray` with shape [y, x, bands]. If your
function has parameters, you can pass in custom values using
a lambda function:
`image_proc_funcs=[lambda x: custom_func(x, param1=10)]`.
show_gdf: geopandas.GeoDataFrame, optional
Vector data (e.g. ESRI shapefiles or GeoJSON) can be optionally
plotted over the top of imagery by supplying a
`geopandas.GeoDataFrame` object. To customise colours used to
plot the vector features, create a new column in the
GeoDataFrame called 'colors' specifying the colour used to plot
each feature: e.g. `gdf['colors'] = 'red'`.
To plot vector features at specific moments in time during the
animation, create new 'start_time' and/or 'end_time' columns in
the GeoDataFrame that define the time range used to plot each
feature. Dates can be provided in any string format that can be
converted using the `pandas.to_datetime()`. e.g.
`gdf['end_time'] = ['2001', '2005-01', '2009-01-01']`
show_date : string or bool, optional
An optional string or bool that defines how (or if) to plot
date annotations for each animation frame. Defaults to
'%d %b %Y'; can be customised to any format understood by
strftime (https://strftime.org/). Set to False to remove date
annotations completely.
show_text : str or list of strings, optional
An optional string or list of strings with a length equal to
the number of timesteps in `ds`. This can be used to display a
static text annotation (using a string), or a dynamic title
(using a list) that displays different text for each timestep.
By default, no text annotation will be plotted.
show_colorbar : bool, optional
An optional boolean indicating whether to include a colourbar
for single-band animations. Defaults to True.
gdf_kwargs : dict, optional
An optional dictionary of keyword arguments to customise the
appearance of a `geopandas.GeoDataFrame` supplied to
`show_gdf`. Keyword arguments are passed to `GeoSeries.plot`
(see http://geopandas.org/reference.html#geopandas.GeoSeries.plot).
For example: `gdf_kwargs = {'linewidth': 2}`.
annotation_kwargs : dict, optional
An optional dict of keyword arguments for controlling the
appearance of text annotations. Keyword arguments are passed
to `matplotlib`'s `plt.annotate`
(see https://matplotlib.org/api/_as_gen/matplotlib.pyplot.annotate.html
for options). For example, `annotation_kwargs={'fontsize':20,
'color':'red', 'family':'serif'}.
imshow_kwargs : dict, optional
An optional dict of keyword arguments for controlling the
appearance of arrays passed to `matplotlib`'s `plt.imshow`
(see https://matplotlib.org/api/_as_gen/matplotlib.pyplot.imshow.html
for options). For example, a green colour scheme and custom
stretch could be specified using:
`onebandplot_kwargs={'cmap':'Greens`, 'vmin':0.2, 'vmax':0.9}`.
(some parameters like 'cmap' will only have an effect for
single-band animations, not three-band RGB animations).
colorbar_kwargs : dict, optional
An optional dict of keyword arguments used to control the
appearance of the colourbar. Keyword arguments are passed to
`matplotlib.pyplot.tick_params`
(see https://matplotlib.org/api/_as_gen/matplotlib.pyplot.tick_params.html
for options). This can be used to customise the colourbar
ticks, e.g. changing tick label colour depending on the
background of the animation:
`colorbar_kwargs={'colors': 'black'}`.
limit: int, optional
An optional integer specifying how many animation frames to
render (e.g. `limit=50` will render the first 50 frames). This
can be useful for quickly testing animations without rendering
the entire time-series.
"""
def _start_end_times(gdf, ds):
"""
Converts 'start_time' and 'end_time' columns in a
`geopandas.GeoDataFrame` to datetime objects to allow vector
features to be plotted at specific moments in time during an
animation, and sets default values based on the first
and last time in `ds` if this information is missing from the
dataset.
"""
# Make copy of gdf so we do not modify original data
gdf = gdf.copy()
# Get min and max times from input dataset
minmax_times = pd.to_datetime(ds.time.isel(time=[0, -1]).values)
# Update both `start_time` and `end_time` columns
for time_col, time_val in zip(['start_time', 'end_time'], minmax_times):
# Add time_col if it does not exist
if time_col not in gdf:
gdf[time_col] = np.nan
# Convert values to datetimes and fill gaps with relevant time value
gdf[time_col] = pd.to_datetime(gdf[time_col], errors='ignore')
gdf[time_col] = gdf[time_col].fillna(time_val)
return gdf
def _add_colorbar(fig, ax, vmin, vmax, imshow_defaults, colorbar_defaults):
"""
Adds a new colorbar axis to the animation with custom minimum
and maximum values and styling.
"""
# Create new axis object for colorbar
cax = fig.add_axes([0.02, 0.02, 0.96, 0.03])
# Initialise color bar using plot min and max values
img = ax.imshow(np.array([[vmin, vmax]]), **imshow_defaults)
fig.colorbar(img,
cax=cax,
orientation='horizontal',
ticks=np.linspace(vmin, vmax, 2))
# Fine-tune appearance of colorbar
cax.xaxis.set_ticks_position('top')
cax.tick_params(axis='x', **colorbar_defaults)
cax.get_xticklabels()[0].set_horizontalalignment('left')
cax.get_xticklabels()[-1].set_horizontalalignment('right')
def _frame_annotation(times, show_date, show_text):
"""
Creates a custom annotation for the top-right of the animation
by converting a `xarray.DataArray` of times into strings, and
combining this with a custom text annotation. Handles cases
where `show_date=False/None`, `show_text=False/None`, or where
`show_text` is a list of strings.
"""
# Test if show_text is supplied as a list
is_sequence = isinstance(show_text, (list, tuple, np.ndarray))
# Raise exception if it is shorter than number of dates
if is_sequence and (len(show_text) == 1):
show_text, is_sequence = show_text[0], False
elif is_sequence and (len(show_text) < len(times)):
raise ValueError(f'Annotations supplied via `show_text` must have '
f'either a length of 1, or a length >= the number '
f'of timesteps in `ds` (n={len(times)})')
times_list = (times.dt.strftime(show_date).values
if show_date else [None] * len(times))
text_list = show_text if is_sequence else [show_text] * len(times)
annotation_list = [
'\n'.join([str(i)
for i in (a, b)
if i])
for a, b in zip(times_list, text_list)
]
return annotation_list
def _update_frames(i, ax, extent, annotation_text, gdf, gdf_defaults,
annotation_defaults, imshow_defaults):
"""
Animation called by `matplotlib.animation.FuncAnimation` to
animate each frame in the animation. Plots array and any text
annotations, as well as a temporal subset of `gdf` data based
on the times specified in 'start_time' and 'end_time' columns.
"""
# Clear previous frame to optimise render speed and plot imagery
ax.clear()
ax.imshow(array[i, ...].clip(0.0, 1.0),
extent=extent,
vmin=0.0,
vmax=1.0,
**imshow_defaults)
# Add annotation text
ax.annotate(annotation_text[i], **annotation_defaults)
# Add geodataframe annotation
if show_gdf is not None:
# Obtain start and end times to filter geodataframe features
time_i = ds.time.isel(time=i).values
# Subset geodataframe using start and end dates
gdf_subset = show_gdf.loc[(show_gdf.start_time <= time_i) &
(show_gdf.end_time >= time_i)]
if len(gdf_subset.index) > 0:
# Set color to geodataframe field if supplied
if ('color' in gdf_subset) and ('color' not in gdf_kwargs):
gdf_defaults.update({'color': gdf_subset['color'].tolist()})
gdf_subset.plot(ax=ax, **gdf_defaults)
# Remove axes to show imagery only
ax.axis('off')
# Update progress bar
progress_bar.update(1)
# Add GeoBox and odc.* accessor to array using `odc-geo`
try:
ds = add_geobox(ds)
except ValueError:
raise ValueError("Unable to determine `ds`'s coordinate "
"reference system (CRS). Please assign a CRS "
"to the array before passing it to this "
"function, e.g.: "
"`ds.odc.assign_crs(crs='EPSG:3577')`")
# Test if bands have been supplied, or convert to list to allow
# iteration if a single band is provided as a string
if bands is None:
raise ValueError(f'Please use the `bands` parameter to supply '
f'a list of one or three bands that exist as '
f'variables in `ds`, e.g. {list(ds.data_vars)}')
elif isinstance(bands, str):
bands = [bands]
# Test if bands exist in dataset
missing_bands = [b for b in bands if b not in ds.data_vars]
if missing_bands:
raise ValueError(f'Band(s) {missing_bands} do not exist as '
f'variables in `ds` {list(ds.data_vars)}')
# Test if time dimension exists in dataset
if 'time' not in ds.dims:
raise ValueError(f"`ds` does not contain a 'time' dimension "
f"required for generating an animation")
# Set default parameters
outline = [PathEffects.withStroke(linewidth=2.5, foreground='black')]
annotation_defaults = {
'xy': (1, 1),
'xycoords': 'axes fraction',
'xytext': (-5, -5),
'textcoords': 'offset points',
'horizontalalignment': 'right',
'verticalalignment': 'top',
'fontsize': 20,
'color': 'white',
'path_effects': outline
}
imshow_defaults = {'cmap': 'magma', 'interpolation': 'nearest'}
colorbar_defaults = {'colors': 'white', 'labelsize': 12, 'length': 0}
gdf_defaults = {'linewidth': 1.5}
# Update defaults with kwargs
annotation_defaults.update(annotation_kwargs)
imshow_defaults.update(imshow_kwargs)
colorbar_defaults.update(colorbar_kwargs)
gdf_defaults.update(gdf_kwargs)
# Get info on dataset dimensions
height, width = ds.odc.geobox.shape
scale = width_pixels / width
left, bottom, right, top = ds.odc.geobox.extent.boundingbox
# Prepare annotations
annotation_list = _frame_annotation(ds.time, show_date, show_text)
# Prepare geodataframe
if show_gdf is not None:
show_gdf = show_gdf.to_crs(ds.odc.geobox.crs)
show_gdf = gpd.clip(show_gdf, mask=box(
left, bottom, right, top)).reindex(show_gdf.index).dropna(how='all')
show_gdf = _start_end_times(show_gdf, ds)
# Convert data to 4D numpy array of shape [time, y, x, bands]
ds = ds[bands].to_array().transpose(..., 'variable')[0:limit, ...]
array = ds.astype(np.float32).values
# Optionally apply image processing along axis 0 (e.g. to each timestep)
bar_format='{l_bar}{bar}| {n_fmt}/{total_fmt} ({remaining_s:.1f} ' \
'seconds remaining at {rate_fmt}{postfix})'
if image_proc_funcs:
print('Applying custom image processing functions')
for i, array_i in tqdm(enumerate(array),
total=len(ds.time),
leave=False,
bar_format=bar_format,
unit=' frames'):
for func in image_proc_funcs:
array_i = func(array_i)
array[i, ...] = array_i
# Clip to percentiles and rescale between 0.0 and 1.0 for plotting
vmin, vmax = np.quantile(array[np.isfinite(array)], q=percentile_stretch)
# Replace with vmin and vmax if present in `imshow_defaults`
if 'vmin' in imshow_defaults:
vmin = imshow_defaults.pop('vmin')
if 'vmax' in imshow_defaults:
vmax = imshow_defaults.pop('vmax')
# Rescale between 0 and 1
array = rescale_intensity(array,
in_range=(vmin, vmax),
out_range=(0.0, 1.0))
array = np.squeeze(array) # remove final axis if only one band
# Set up figure
fig, ax = plt.subplots()
fig.set_size_inches(width * scale / 72, height * scale / 72, forward=True)
fig.subplots_adjust(left=0, bottom=0, right=1, top=1, wspace=0, hspace=0)
# Optionally add colorbar
if show_colorbar & (len(bands) == 1):
_add_colorbar(fig, ax, vmin, vmax, imshow_defaults, colorbar_defaults)
# Animate
print(f'Exporting animation to {output_path}')
anim = FuncAnimation(
fig=fig,
func=_update_frames,
fargs=(
ax, # axis to plot into
[left, right, bottom, top], # imshow extent
annotation_list, # list of text annotations
show_gdf, # geodataframe to plot over imagery
gdf_defaults, # any kwargs used to plot gdf
annotation_defaults, # kwargs for annotations
imshow_defaults), # kwargs for imshow
frames=len(ds.time),
interval=interval,
repeat=False)
# Set up progress bar
progress_bar = tqdm(total=len(ds.time),
unit=' frames',
bar_format=bar_format)
# Export animation to file
if Path(output_path).suffix == '.gif':
anim.save(output_path, writer='pillow')
else:
anim.save(output_path, dpi=72)
# Update progress bar to fix progress bar moving past end
if progress_bar.n != len(ds.time):
progress_bar.n = len(ds.time)
progress_bar.last_print_n = len(ds.time)
def _degree_to_zoom_level(l1, l2, margin=0.0):
"""
Helper function to set zoom level for `display_map`
"""
degree = abs(l1 - l2) * (1 + margin)
zoom_level_int = 0
if degree != 0:
zoom_level_float = math.log(360 / degree) / math.log(2)
zoom_level_int = int(zoom_level_float)
else:
zoom_level_int = 18
return zoom_level_int
[docs]
def plot_wo(wo, legend=True, **plot_kwargs):
"""Plot a water observation bit flag image.
Parameters
----------
wo : xr.DataArray
A DataArray containing water observation bit flags.
legend : bool
Whether to plot a legend. Default True.
plot_kwargs : dict
Keyword arguments passed on to DataArray.plot.
Returns
-------
plot
"""
cmap = mcolours.ListedColormap([
np.array([150, 150, 110]) / 255, # dry - 0
np.array([0, 0, 0]) / 255, # nodata, - 1
np.array([119, 104, 87]) / 255, # terrain - 16
np.array([89, 88, 86]) / 255, # cloud_shadow - 32
np.array([216, 215, 214]) / 255, # cloud - 64
np.array([242, 220, 180]) / 255, # cloudy terrain - 80
np.array([79, 129, 189]) / 255, # water - 128
np.array([51, 82, 119]) / 255, # shady water - 160
np.array([186, 211, 242]) / 255, # cloudy water - 192
])
bounds = [
0,
1,
16,
32,
64,
80,
128,
160,
192,
255,
]
norm = mcolours.BoundaryNorm(np.array(bounds) - 0.1, cmap.N)
cblabels = [
'dry', 'nodata', 'terrain', 'cloud shadow', 'cloud', 'cloudy terrain',
'water', 'shady water', 'cloudy water'
]
try:
im = wo.plot.imshow(cmap=cmap,
norm=norm,
add_colorbar=legend,
**plot_kwargs)
except AttributeError:
im = wo.plot(cmap=cmap, norm=norm, add_colorbar=legend, **plot_kwargs)
if legend:
try:
cb = im.colorbar
except AttributeError:
cb = im.cbar
ticks = cb.get_ticks()
cb.set_ticks(ticks[:-1] + np.diff(ticks) / 2)
cb.set_ticklabels(cblabels)
return im
[docs]
def plot_fmask(fmask, legend=True, **plot_kwargs):
"""
Plot an enumerated FMask flag image with human-readable colours.
Parameters
----------
fmask : xr.DataArray
A DataArray containing Fmask flags.
legend : bool
Whether to plot a legend. Default True.
plot_kwargs : dict
Keyword arguments passed on to DataArray.plot.
Returns
-------
plot
"""
cmap = mcolours.ListedColormap([
np.array([0, 0, 0]) / 255, # nodata - 0
np.array([132, 162, 120]) / 255, # clear - 1
np.array([208, 207, 206]) / 255, # cloud - 2
np.array([70, 70, 51]) / 255, # cloud_shadow - 3
np.array([224, 237, 255]) / 255, # snow - 4
np.array([71, 91, 116]) / 255, # water - 5
])
bounds = [0, 1, 2, 3, 4, 5, 6]
norm = mcolours.BoundaryNorm(np.array(bounds) - 0.1, cmap.N)
cblabels = ['nodata', 'clear', 'cloud', 'shadow', 'snow', 'water']
try:
im = fmask.plot.imshow(cmap=cmap,
norm=norm,
add_colorbar=legend,
**plot_kwargs)
except AttributeError:
im = fmask.plot(cmap=cmap,
norm=norm,
add_colorbar=legend,
**plot_kwargs)
if legend:
try:
cb = im.colorbar
except AttributeError:
cb = im.cbar
ticks = cb.get_ticks()
cb.set_ticks(ticks[:-1] + np.diff(ticks) / 2)
cb.set_ticklabels(cblabels)
return im
[docs]
def plot_variable_images(img_collection):
"""
Plot a dynamic number of images from a xarray dataset that
includes Date and Index in the title. Optional ability to
also include the sensor in the title if a 'sensor' attribute
is added to the dataset using dataset.assign_attrs
Parameters
----------
img_collection : xr.Dataset
A Dataset containing imagery with RBG bands
Returns
-------
plot
"""
# Check that img_collection is a xarray dataset
if not isinstance(img_collection, xr.Dataset):
raise TypeError("img_collection must be a xarray dataset.")
# Calculate number of images in `img_collection`
plot_count = img_collection.dims["time"]
# Check if dataset has 0 images
if plot_count == 0:
if hasattr(img_collection, "sensor"):
raise ValueError("The {} dataset has no images to display for the "
"given query parameters".format(img_collection.sensor))
else:
raise ValueError("The supplied xarray dataset has no images to "
"display for the given query parameters")
# Divide the number of images by 2 rounding up to calculate the
# number of rows for the below figure are needed
plot_rows = math.ceil(plot_count / 2)
# Construct a figure to visualise the imagery
f, axarr = plt.subplots(plot_rows, 2, figsize=(10, plot_rows * 4.5),
squeeze=False)
# Flatten the subplots so they can easily be enumerated through
axarr = axarr.flatten()
# Iterate through each image in the dataset and plot
# each image as a RGB on a subplot
for t in range(plot_count):
rgb(
img_collection.isel(time=t),
bands=["nbart_red", "nbart_green", "nbart_blue"],
ax=axarr[t],
robust=True,
)
# Test to see if the dataset has a custom 'sensor' attribute.
# If so, include the string in each subplot title.
if hasattr(img_collection, "sensor"):
title = (
str(img_collection.time[t].values)[:10]
+ " || Index: "
+ str(t)
+ " || Sensor: "
+ img_collection.sensor
)
else:
title = (
str(img_collection.time[t].values)[:10]
+ " || Index: "
+ str(t)
)
# Set subplot title, axis label, and shrink offset text
axarr[t].set_title(title)
axarr[t].set_xlabel("X coordinate")
axarr[t].set_ylabel("Y coordinate")
axarr[t].yaxis.offsetText.set_fontsize(6)
axarr[t].xaxis.offsetText.set_fontsize(6)
# Adjust padding arround subplots to prevent overlapping elements
plt.tight_layout()
# Remove the last subplot if an odd number of images are being displayed
if plot_count % 2 != 0:
f.delaxes(axarr[plot_count])
