DEA Land Cover pixel drill
Sign up to the DEA Sandbox to run this notebook interactively from a browser
Compatibility: Notebook currently compatible with the
DEA Sandbox
environmentProducts used: ga_ls_landcover_class_cyear_2
Background
Land cover is the physical surface of the Earth, including trees, shrubs, grasses, soils, exposed rocks, water bodies, plantations, crops and built structures. including trees, shrubs, grasses, soils, exposed rocks, water bodies, plantations, crops and built structures.
The Digital Earth Australia (DEA) Land Cover product classifies the Australian continent’s surface using the FAO LCCS taxonomy, providing detailed information on the presence of water, vegetation cover and the urban environment. For a more in-depth introduction to the DEA Land Cover dataset, see the DEA Land Cover notebook and the DEA Land Cover product details page.
Trends in land cover over time provide important information about environmental conditions and human activity. For example, measuring change in land cover over time can describe the natural cycles of greening and drying in an ephemeral wetland, or show how an ecosystem responds after a bushfire.
Description
A pixel drill is a method of analysing time-deep satellite data; a single spatial point is chosen and the characteristics of that point through time are plotted. This notebook conducts a pixel drill on the Level 3 classification layer of the DEA Land Cover product, visualising how this location on Australia’s surface has or has not changed over 30+ years.
Select an area of interest and load data.
Visualise the data and select point for pixel drill.
Conduct pixel drill.
Plot land cover history of the chosen point.
Getting started
To run this analysis, run all the cells in the notebook, starting with the ‘Load packages and connect to the datacube’ cell.
Load packages
[1]:
%matplotlib widget
import os
import sys
import datacube
import matplotlib.pyplot as plt
sys.path.insert(1, os.path.abspath("../Tools"))
from dea_tools.landcover import lc_colourmap
from dea_tools.plotting import display_map
Connect to the datacube
Connect to the datacube so we can access DEA data.
[2]:
dc = datacube.Datacube(app="Land_cover_pixel_drill")
Define custom functions
This notebook uses the following custom functions to format the data and create the interactive plot.
[3]:
def rescale_level3(data):
"""
Rescale distribution of Level 3 classes to be 0-6 for neat plotting in pixel drill plot.
Input: xarray DataSet or DataArray containing Level 3 land cover data.
"""
scale_and_lables = {0: 0, 111: 1, 112: 2, 124: 3, 215: 4, 216: 5, 220: 6}
for key, value in scale_and_lables.items():
data = data.where(data != key, value)
return data
def click_for_coords(event):
"""
Define behaviour widget to enable extraction of x,y values for point manually selected with mouse click on interactive plot.
For use with Matplotlib Imshow.
The x and y returned are the coordinate positions of the pixel in the x and y axis, rather than latitude/longitude values of the point.
The coordinate positions of the pixel this is what you need later to conduct a pixel drill.
Returns: ix, iy = position of selected pixel in the x and y coordinates of the array.
"""
global ix, iy
ix, iy = event.xdata, event.ydata
global coords
coords.append((ix, iy))
fig.canvas.mpl_disconnect(cid)
return
Define area of interest and load data
[4]:
# Select central coordinates for an area of interest
# Lake George
point_x, point_y = (-35.1064, 149.4248)
# Combine central coordinates with buffer values to create the latitude and longitude range
lat = (point_x - 0.1, point_x + 0.1)
lon = (point_y - 0.1, point_y + 0.1)
# Display area on map
display_map(x=lon, y=lat)
[4]:
[5]:
# Create the 'query' dictionary object, which contains the coordinates defined above
query = {
"y": lat,
"x": lon,
"time": ("1988", "2020"),
}
# Load DEA Land Cover data from the datacube
all_data = dc.load(
product="ga_ls_landcover_class_cyear_2",
output_crs="EPSG:3577",
measurements=["level3"],
resolution=(-25, 25),
**query
)
[6]:
# Have a look at the data we loaded
all_data
[6]:
<xarray.Dataset> Dimensions: (time: 33, y: 982, x: 842) Coordinates: * time (time) datetime64[ns] 1988-01-01 1989-01-01 ... 2020-01-01 * y (y) float64 -3.93e+06 -3.93e+06 ... -3.954e+06 -3.954e+06 * x (x) float64 1.569e+06 1.569e+06 1.569e+06 ... 1.59e+06 1.59e+06 spatial_ref int32 3577 Data variables: level3 (time, y, x) uint8 112 112 112 112 111 ... 112 112 112 112 112 Attributes: crs: EPSG:3577 grid_mapping: spatial_ref
- time: 33
- y: 982
- x: 842
- time(time)datetime64[ns]1988-01-01 ... 2020-01-01
- units :
- seconds since 1970-01-01 00:00:00
array(['1988-01-01T00:00:00.000000000', '1989-01-01T00:00:00.000000000', '1990-01-01T00:00:00.000000000', '1991-01-01T00:00:00.000000000', '1992-01-01T00:00:00.000000000', '1993-01-01T00:00:00.000000000', '1994-01-01T00:00:00.000000000', '1995-01-01T00:00:00.000000000', '1996-01-01T00:00:00.000000000', '1997-01-01T00:00:00.000000000', '1998-01-01T00:00:00.000000000', '1999-01-01T00:00:00.000000000', '2000-01-01T00:00:00.000000000', '2001-01-01T00:00:00.000000000', '2002-01-01T00:00:00.000000000', '2003-01-01T00:00:00.000000000', '2004-01-01T00:00:00.000000000', '2005-01-01T00:00:00.000000000', '2006-01-01T00:00:00.000000000', '2007-01-01T00:00:00.000000000', '2008-01-01T00:00:00.000000000', '2009-01-01T00:00:00.000000000', '2010-01-01T00:00:00.000000000', '2011-01-01T00:00:00.000000000', '2012-01-01T00:00:00.000000000', '2013-01-01T00:00:00.000000000', '2014-01-01T00:00:00.000000000', '2015-01-01T00:00:00.000000000', '2016-01-01T00:00:00.000000000', '2017-01-01T00:00:00.000000000', '2018-01-01T00:00:00.000000000', '2019-01-01T00:00:00.000000000', '2020-01-01T00:00:00.000000000'], dtype='datetime64[ns]')
- y(y)float64-3.93e+06 -3.93e+06 ... -3.954e+06
- units :
- metre
- resolution :
- -25.0
- crs :
- EPSG:3577
array([-3929587.5, -3929612.5, -3929637.5, ..., -3954062.5, -3954087.5, -3954112.5])
- x(x)float641.569e+06 1.569e+06 ... 1.59e+06
- units :
- metre
- resolution :
- 25.0
- crs :
- EPSG:3577
array([1569212.5, 1569237.5, 1569262.5, ..., 1590187.5, 1590212.5, 1590237.5])
- spatial_ref()int323577
- spatial_ref :
- PROJCS["GDA94 / Australian Albers",GEOGCS["GDA94",DATUM["Geocentric_Datum_of_Australia_1994",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],AUTHORITY["EPSG","6283"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4283"]],PROJECTION["Albers_Conic_Equal_Area"],PARAMETER["latitude_of_center",0],PARAMETER["longitude_of_center",132],PARAMETER["standard_parallel_1",-18],PARAMETER["standard_parallel_2",-36],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","3577"]]
- grid_mapping_name :
- albers_conical_equal_area
array(3577, dtype=int32)
- level3(time, y, x)uint8112 112 112 112 ... 112 112 112 112
- units :
- 1
- nodata :
- 0
- crs :
- EPSG:3577
- grid_mapping :
- spatial_ref
array([[[112, 112, 112, ..., 112, 111, 111], [112, 112, 112, ..., 112, 112, 112], [112, 112, 112, ..., 111, 111, 111], ..., [112, 112, 112, ..., 112, 112, 112], [112, 112, 112, ..., 112, 112, 112], [112, 112, 112, ..., 112, 112, 112]], [[112, 112, 112, ..., 112, 112, 112], [112, 112, 112, ..., 111, 112, 111], [112, 112, 112, ..., 111, 111, 111], ..., [112, 112, 112, ..., 112, 112, 112], [112, 112, 112, ..., 112, 112, 112], [112, 112, 112, ..., 112, 112, 112]], [[112, 112, 112, ..., 111, 112, 111], [112, 112, 112, ..., 111, 112, 111], [112, 112, 112, ..., 111, 111, 111], ..., ... ..., [112, 112, 112, ..., 112, 112, 112], [112, 112, 112, ..., 112, 112, 112], [112, 112, 112, ..., 112, 112, 112]], [[112, 112, 112, ..., 112, 112, 112], [112, 112, 112, ..., 112, 112, 112], [112, 112, 112, ..., 112, 112, 112], ..., [112, 112, 112, ..., 112, 112, 112], [112, 112, 112, ..., 112, 112, 112], [112, 112, 112, ..., 112, 112, 112]], [[112, 112, 112, ..., 111, 112, 111], [112, 112, 112, ..., 111, 111, 111], [112, 112, 112, ..., 111, 111, 111], ..., [112, 112, 112, ..., 112, 112, 112], [112, 112, 112, ..., 112, 112, 112], [112, 112, 112, ..., 112, 112, 112]]], dtype=uint8)
- crs :
- EPSG:3577
- grid_mapping :
- spatial_ref
[7]:
# Select Level 3
level3 = all_data.level3
Visualise the data and select point for pixel drill
This plot has an interactive function enabled which allows users to select a specific point to analyse in greater detail further down the notebook. To use this feature, select a point on the plot you wish to interrogate by clicking on it with your mouse. Then run the next code cell below. The point you selected will then be highlighted on this plot with a red dot. The most recent year has been set as default to visualise in this step.
You can move the dot by clicking on a new location and then re-running the next cell after the interactive plot again. Unless you have since completed the notebook and generated the graph at the end, in which case you need to re-run this code first in order to re-generate the interactive plot before selecting a new point. If you click multiple times on the plot before running the next piece of code only the most recent location will be used. If you run into trouble, restarting the kernel of the notebook and running from the first cell again may help.
[ ]:
# Build the colour map for level 3
layer_cmap, layer_norm = lc_colourmap("level3")
# Select year to plot
plot_year = "2020"
# Create interactive plot
fig, ax = plt.subplots(figsize=(6, 6))
ax = plt.imshow(
level3.sel(time=plot_year, method="nearest").squeeze(),
cmap=layer_cmap,
norm=layer_norm,
interpolation="nearest",
)
cid = fig.canvas.mpl_connect("button_press_event", click_for_coords)
plt.title("Interactive Plot")
plt.show()
Run next cell to enter or update selected point
[9]:
# Enter selected point
plt.clf()
ax = plt.imshow(
level3.sel(time=plot_year, method="nearest").squeeze(),
cmap=layer_cmap,
norm=layer_norm,
interpolation="nearest",
)
plt.scatter(x=ix, y=iy, c="r", s=80)
[9]:
<matplotlib.collections.PathCollection at 0x7fbd51522fd0>
Troubleshooting:
If an error is appearing above this test saying NameError: name ‘ix’ is not defined you have not chosen a point to conduct the pixel drill on or entered it. Please select a point on the interactive plot and run the cell above to enter the point.
Conduct pixel drill
[10]:
# Find nearest whole number value for x,y
xcoord = int(ix)
ycoord = int(iy)
# Conduct pixel drill on loaded dataset
pixel_timeseries = level3.isel(
x=xcoord,
y=ycoord,
)
# Rescale the range of Level 3 classes for beautiful plotting
new_data = rescale_level3(pixel_timeseries)
Plot land cover history of the chosen point
To visualise the results of the pixel drill we will generate a line graph with the time on the x axis and the land cover classes on the y axis. The following code also re-labels the y axis with the proper names of the land cover classes to make the plot easily readable.
[ ]:
# Define new figure and plot
fig, ax = plt.subplots(figsize=(8, 4))
# Relabel y axis with Level 3 class names
number_labels = (0, 1, 2, 3, 4, 5, 6)
label_classes = (
"No Data",
"Cultivated Terrestrial\n Vegetation",
"Natural Terrestrial\n Vegetation",
"Natural Aquatic\n Vegetation",
"Artificial Surface",
"Natural Bare\n Surface",
"Water",
)
plt.yticks(number_labels, label_classes)
# Add data to plot
new_data.plot(linestyle="--", c="b", marker="8", mec="b", mfc="b")
plt.xlabel("year")
plt.ylabel("Level 3 Land Cover class")
plt.title("Land Cover History for Selected Point")
plt.grid()
plt.tight_layout()
plt.show
Additional information
License: The code in this notebook is licensed under the Apache License, Version 2.0. Digital Earth Australia data is licensed under the Creative Commons by Attribution 4.0 license.
Contact: If you need assistance, please post a question on the Open Data Cube Discord chat or on the GIS Stack Exchange using the open-data-cube
tag (you can view previously asked questions here). If you would like to report an issue with this notebook, you can file one on
GitHub.
Last modified: December 2023
Compatible datacube version:
[12]:
print(datacube.__version__)
1.8.6