Maintainer: Xavier Ho xavier.ho@csiro.au
Please contact Xavier for any queries about this example.
This example uses the Analytics Execution Engine of AGDC v2 to query Landsat imagery. Two bands are combined to calculate an average, and the result is saved to an image.
In [3]:
from pprint import pprint
from datetime import datetime
import xarray as xr
import matplotlib
import matplotlib.image
%matplotlib inline
import datacube
from datacube.api import API, geo_xarray
from datacube.analytics.analytics_engine import AnalyticsEngine
from datacube.execution.execution_engine import ExecutionEngine
from datacube.analytics.utils.analytics_utils import plot
print('This example runs on Data Cube v2/{}.'.format(datacube.__version__))
First, we make a query to the datacube to find out what datasets we have.
In [4]:
dc_a = AnalyticsEngine()
dc_e = ExecutionEngine()
dc_api = API()
print(dc_api.list_field_values('product')) # 'LEDAPS' should be in the list
print(dc_api.list_field_values('platform')) # 'LANDSAT_5' should be in the list
Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) is a NASA-funded project to map North American forest disturbance since 1975. We have datasets in the same format for Australia.
Let's find out what kind of datasets we have for Landsat 5.
In [6]:
query = {
'product': 'LEDAPS',
'platform': 'LANDSAT_5',
}
descriptor = dc_api.get_descriptor(query, include_storage_units=False)
pprint(descriptor)
For Landsat 5, Band 1-3 are Blue, Green, and Red visible light spectrum bands.
To set up the Engine, we first need to instantiate the modules and setup query parameters.
create_arraysets up the platform and product we are interested in querying, as well as the bands (variables) of the satellite data set. We also limit the amount of data processed by a long-lat boundary and time.
apply_expressionbinds the variables into a generic string to execute, in this case an average of two bands.
execute_planis when the computation is actually run and returned.
In [13]:
dimensions = {
'x': {
'range': (140, 141)
},
'y': {
'range': (-35.5, -36.5)
},
'time': {
'range': (datetime(2011, 10, 17), datetime(2011, 10, 18))
}
}
red = dc_a.create_array(('LANDSAT_5', 'LEDAPS'), ['band3'], dimensions, 'red')
green = dc_a.create_array(('LANDSAT_5', 'LEDAPS'), ['band2'], dimensions, 'green')
blue = dc_a.create_array(('LANDSAT_5', 'LEDAPS'), ['band1'], dimensions, 'blue')
Now we have created references to the green and blue bands, we can do simple band maths.
In [19]:
blue_result = dc_a.apply_expression([blue], 'array1', 'blue')
dc_e.execute_plan(dc_a.plan)
plot(dc_e.cache['blue'])
In [22]:
turbidity = dc_a.apply_expression([blue, green, red], '(array1 + array2 - array3) / 2', 'turbidity')
dc_e.execute_plan(dc_a.plan)
plot(dc_e.cache['turbidity'])
geo_xarray.reproject reprojects northings and eastings to longitude and latitude units.
Let's reproject the file into the common longitude-latitude projection, and save it to a picture.
In [23]:
result = dc_e.cache['turbidity']['array_result']['turbidity']
reprojected = datacube.api.geo_xarray.reproject(result.isel(time=0), 'EPSG:3577', 'WGS84')
pprint(reprojected)
reprojected.plot.imshow()
matplotlib.image.imsave('turbidity.png', reprojected)
The boundaries in long-lat are as follows:
In [24]:
map(float, (reprojected.x[0], reprojected.x[-1], reprojected.y[0], reprojected.y[-1]))
Out[24]:
In [ ]: