<jupyter_start><jupyter_text>Introduction to plotting In this notebook we cover `xarray` data plotting functionality for generating simple and quick plots from multi-dimensional datasets.  The DEA Python library also provides a custom plotting module with additional convenient functionality.<jupyter_code>%matplotlib inline

import datacube
from check_answer import check_answer
import sys

sys.path.append("./Scripts")
from dea_plotting import display_map
from dea_plotting import rgb

dc = datacube.Datacube(app="Plotting")<jupyter_output><empty_output><jupyter_text>The `display_map()` function is provided in the DEA environment for interactively visualising the geographical extents of a query. It uses the Leaflet web mapping library presenting the output inside Jupyter notebooks. This function provides a convenient way of checking that the boundaries in your query correspond to the right location.<jupyter_code>lat_range = (-27.58, -27.666)
lon_range = (153.3, 153.4)

display_map(x=lon_range, y=lat_range)<jupyter_output><empty_output><jupyter_text>The variables specified above can used in a query to load data from the DEA datacube using `dc.load()`.<jupyter_code>ds = dc.load(product="ls8_nbart_geomedian_annual",
             x=lon_range,
             y=lat_range,
             time=("2015", "2019"))

ds
da = ds.swir1.sel(time='2017')

da<jupyter_output><empty_output><jupyter_text>Exercise 3.1: If you remember from the XArray tutorial, there are two different ways of selecting parts of Datasets and DataArrays. First uses `.sel()` to when we want to pass values for the dimensions and `.isel()` when we want to pass positional indexes. Can you request the same temporal slice as in the example above but using `.isel()` this time?<jupyter_code>answ = ds.?.isel(?)

check_answer("DEA3.1", str(answ.time.values)[:4])<jupyter_output><empty_output><jupyter_text>The plotting functionality that we saw in the XArray tutorial also applies to DEA Datasets:>Tip: Remember that for plotting we always need to select one temporal element from the Dataset.<jupyter_code>ds.swir1.isel(time=0).plot()<jupyter_output><empty_output><jupyter_text>Exercise 3.2: Can you create a plot for the red channel in 2018?  It is often useful to produce plots for a single measurement across time, for example to compare change between satellite observations or summary datasets. To plot multiple images, skip the `isel()` step above and plot the entire `xarray.DataArray` directly.  To plot multiple timesteps in one figure, it is necessary to instruct the `.plot()` function to put each timestep in a different column. This is done by specifying `.plot(col="time")`:<jupyter_code>ds.swir1.plot(col="time")<jupyter_output><empty_output><jupyter_text>>Note: This kind of plotting is called "facetted plotting". For more information, refer to the [xarray documentation](http://xarray.pydata.org/en/stable/plotting.htmlfaceting)  The plots above are dark and difficult to see details in them. To improve the appearance of `xarray` plots, use the `robust=True` argument to optimise the plot colours by clipping extreme values or outliers. This will use the 2nd and 98th percentiles of the data to compute the color limits:<jupyter_code>ds.swir1.plot(col="time", robust=True)<jupyter_output><empty_output><jupyter_text>Plots can be further customised by adding custom colour maps/styles using the `cmap` parameter. When choosing a colour map for a plot, it is important to choose a set of colours that are perceived logically by the human eye. The best colour maps are "perceptually uniform": these colour maps increase logically from dark to light colours, where equal increases in lightness/darkness correspond to equal changes in data values. Some best-practice perceptually uniform colour maps include:```"viridis", "plasma", "inferno", "magma", "cividis"```> For further reading about perceptually uniform colour maps in data visualisation, refer to the [matplotlib documentation](https://matplotlib.org/3.1.0/tutorials/colors/colormaps.html) It is also important to consider colour blindness when selecting a colour map. `XArray` supports many colour maps from the "colorbrewer" family of colour maps which are optimised for colour blindness. You can use the interactive [online tool](http://colorbrewer2.org) to browse all available colour maps, or choose from one of the following commonly used options:```"Greys", "Purples", "Blues", "Greens", "Oranges", "Reds","YlOrBr", "YlOrRd", "OrRd", "PuRd", "RdPu", "BuPu", "GnBu", "PuBu", "YlGnBu", "PuBuGn", "BuGn", "YlGn"``` For a full list of available colour maps you can refer to [this list](https://matplotlib.org/3.1.0/tutorials/colors/colormaps.html). The example cell below plots the data with the perceptually uniform `magma` colour map:<jupyter_code>ds.swir1.plot(col="time", robust=True, cmap="magma")<jupyter_output><empty_output><jupyter_text>Exercise 3.3: Experiment with the previous colour maps and create different displays with the contents of the dataset.  Although `XArray` makes it easy to plot single band images, plotting a three band RGB colour image is less straightforward. If you recall the not easy to remember chain of methods that we have been using previously:<jupyter_code>ds[['red', 'green', 'blue']].isel(time=0).to_array().plot.imshow(robust=True, figsize=(6,6))<jupyter_output><empty_output><jupyter_text>To make the previous process easier, the [dea-notebooks repository](https://github.com/GeoscienceAustralia/dea-notebooks) provides a custom `rgb()` function that is designed for plotting three band images. The `rgb()` function maps three data variables/measurements from the loaded dataset to the red, green and blue channels that are used to make a three-colour image. The `rgb()` It requires the minimum input of:* `ds:` The `xarray.Dataset` object* `bands:` Three bands for display (these must be measurements found in the dataset) * `index:` The timestep to view, default is `0` For example for plotting a single timestep:<jupyter_code># View a red, green, blue (true colour) image of the first timestep
rgb(ds, bands=["red", "green", "blue"], index=0)<jupyter_output><empty_output><jupyter_text>Exercise 3.4: Can you create a false colour image, with the following mapping (`swir1`, `nir`, `green`)?<jupyter_code>answ = [?,?,?]

rgb(ds, bands=answ, index=0)

check_answer("DEA4.4", answ)<jupyter_output><empty_output><jupyter_text>The `rgb()` function can also do display multiple timesteps, as long as a list of timesteps to view is provided to the `index` argument, e.g. `index=[X1, X2, ...]`. The example cell below plots the first and fifth image in the dataset using `index=[0, 4]` (remembering that counting in Python starts at 0):<jupyter_code>rgb(ds, bands=['red', 'green', 'blue'], index=[0, 3])<jupyter_output><empty_output><jupyter_text>It is also possible to use `rgb()` to plot all timesteps in a dataset using the `col="time"` syntax demonstrated in the [single band example above](Plotting-multiple-timesteps):<jupyter_code>rgb(ds, bands=['red', 'green', 'blue'], col="time")<jupyter_output><empty_output><jupyter_text>By default, `rgb()` generates plots with `robust=True` to improve the appearance of the images by clipping out the darkest and brightest 2% of pixels, using the 2nd and 98th percentiles of the data to compute the colour limits If this default provides poor results, the plot's colour stretch can be customised using the `percentile_stretch` parameter. This clips the most extreme minimum and maximum values in the dataset, improving the contrast and appearance of the plot. For example, specifying `percentile_stretch=[0.05, 0.95]` will clip out the darkest and brightest 5% of pixels, focusing the colour stretch on the remaining 90% of less extreme values:<jupyter_code>rgb(ds, 
    bands=['red', 'green', 'blue'], 
    index=0, 
    percentile_stretch=[0.05, 0.95])<jupyter_output><empty_output><jupyter_text>Exercise 3.5: Can you create an NDVI composite from the previous dataset and use the `.plot()` using the `summer_r` cmap and passing `'time'` to the col parameter?<jupyter_code>ndvi = ?

ndvi.plot(?, ?)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09676629103380696, 'fraction_numerical': 0.01212640862322391, 'mean_word_length': 5.919491525423729, 'pattern_counts': {'":': 1, '<': 61, '<?xml version=': 0, '>': 64, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '115364', 'license': 'Apache-2.0', 'n_tokens_mistral': 2493, 'n_tokens_neox': 2360, 'n_words': 968, 'path': 'dea_materials/day2/3_dea_plotting.ipynb', 'repo_name': 'ANU-WALD/dea_course', 'size': '8164'}