Working with time in xarray¶
Products used: s2_l2a
Keywords analysis; time series, data used; sentinel-2, data methods; groupby,:index:data methods; nearest, index:data methods; interpolating, data methods; resampling, data methods; compositing
Background¶
Time series data is a series of data points usually captured at successively spaced points in time. In a remote-sensing context, time series data is a sequence of discrete satellite images taken at the same area at successive times. Time series analysis uses different methods to extract meaningful statistics, patterns and other characteristics of the data. Time series data and analysis has widespread application ranging from monitoring agricultural crops, natural vegetation change detection, mineral prospectivity mapping, and tidal height modelling.
Description¶
The xarray
Python package provides many useful techniques for dealing with time series data that can be applied to Digital Earth Africa data. This notebook demonstrates how to use xarray
techniques to:
Select different time periods of data (e.g. year, month, day) from an
xarray.Dataset
Use datetime accessors to extract additional information from a dataset’s
time
dimensionSummarise time series data for different time periods using
.groupby()
and.resample()
Interpolate time series data to estimate landscape conditions at a specific date that the satellite did not observe
For additional information about the techniques demonstrated below, refer to the xarray
time series data guide.
Getting started¶
To run this analysis, run all the cells in the notebook, starting with the “Load packages” cell.
Load packages¶
[1]:
%matplotlib inline
import datacube
import matplotlib.pyplot as plt
import numpy as np
from deafrica_tools.datahandling import load_ard, mostcommon_crs
Connect to the datacube¶
[2]:
dc = datacube.Datacube(app='Working_with_time')
Loading Landsat data¶
First, we load in around two years’ of Sentinel-2 data, using the load_ard function and filtering for timesteps with at least 95% good-quality pixels.
[3]:
lat, lon = 13.94, -16.54
buffer = 0.125
# Create a reusable query
query = {
'x': (lon-buffer, lon+buffer),
'y': (lat+buffer, lat-buffer),
'time': ('2018-01', '2019-12'),
'resolution': (-20, 20),
'measurements':['red', 'green', 'blue', 'nir']
}
# Identify the most common projection system in the input query
output_crs = mostcommon_crs(dc=dc, product='s2_l2a', query=query)
# Load available data from Landsat 8 and filter to retain only times
# with at least 95% good data
ds = load_ard(dc=dc,
products=['s2_l2a'],
min_gooddata=0.95,
output_crs=output_crs,
align=(15, 15),
**query)
Using pixel quality parameters for Sentinel 2
Finding datasets
s2_l2a
Counting good quality pixels for each time step
Filtering to 34 out of 144 time steps with at least 95.0% good quality pixels
Applying pixel quality/cloud mask
Loading 34 time steps
Explore xarray data using time¶
Here we will explore several ways to utilise the time dimension within an xarray.Dataset
. This section outlines selecting, summarising and interpolating data at specific times.
Indexing by time¶
We can select data for an entire year by passing a string to .sel()
:
[4]:
ds.sel(time='2018')
[4]:
<xarray.Dataset> Dimensions: (time: 12, x: 1361, y: 1392) Coordinates: * time (time) datetime64[ns] 2018-01-08T11:46:56 ... 2018-12-14T11:... * y (y) float64 1.556e+06 1.556e+06 ... 1.528e+06 1.528e+06 * x (x) float64 3.2e+05 3.2e+05 3.201e+05 ... 3.472e+05 3.472e+05 spatial_ref int32 32628 Data variables: red (time, y, x) float32 2254.0 1988.0 2706.0 ... 1732.0 1747.0 green (time, y, x) float32 1502.0 1348.0 1825.0 ... 1328.0 1336.0 blue (time, y, x) float32 815.0 742.0 1010.0 ... 1016.0 1002.0 nir (time, y, x) float32 3789.0 3519.0 3973.0 ... 2483.0 2496.0 Attributes: crs: epsg:32628 grid_mapping: spatial_ref
- time: 12
- x: 1361
- y: 1392
- time(time)datetime64[ns]2018-01-08T11:46:56 ... 2018-12-...
- units :
- seconds since 1970-01-01 00:00:00
array(['2018-01-08T11:46:56.000000000', '2018-01-23T11:41:48.000000000', '2018-02-17T11:37:13.000000000', '2018-02-22T11:42:49.000000000', '2018-03-04T11:47:19.000000000', '2018-03-09T11:46:52.000000000', '2018-04-03T11:47:15.000000000', '2018-05-23T11:39:00.000000000', '2018-10-15T11:35:46.000000000', '2018-11-19T11:47:30.000000000', '2018-12-09T11:47:28.000000000', '2018-12-14T11:47:31.000000000'], dtype='datetime64[ns]')
- y(y)float641.556e+06 1.556e+06 ... 1.528e+06
- units :
- metre
- resolution :
- -20.0
- crs :
- epsg:32628
array([1555545., 1555525., 1555505., ..., 1527765., 1527745., 1527725.])
- x(x)float643.2e+05 3.2e+05 ... 3.472e+05
- units :
- metre
- resolution :
- 20.0
- crs :
- epsg:32628
array([320025., 320045., 320065., ..., 347185., 347205., 347225.])
- spatial_ref()int3232628
- spatial_ref :
- PROJCS["WGS 84 / UTM zone 28N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-15],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32628"]]
- grid_mapping_name :
- transverse_mercator
array(32628, dtype=int32)
- red(time, y, x)float322254.0 1988.0 ... 1732.0 1747.0
- units :
- 1
- crs :
- epsg:32628
- grid_mapping :
- spatial_ref
array([[[2254., 1988., 2706., ..., 2111., 2346., 1866.], [2002., 2137., 2540., ..., 1565., 1922., 1982.], [1552., 1970., 2724., ..., 1982., 1975., 1932.], ..., [ 377., 386., 387., ..., 1710., 1561., 1535.], [ 367., 373., 378., ..., 1568., 1507., 1581.], [ 392., 391., 372., ..., 1597., 1602., 1586.]], [[2335., 2147., 2643., ..., 2094., 2332., 2069.], [2016., 2208., 2517., ..., 1813., 2144., 2168.], [1783., 2293., 2686., ..., 2256., 2145., 2113.], ..., [ 578., 588., 592., ..., 1807., 1723., 1729.], [ 584., 593., 585., ..., 1762., 1710., 1738.], [ 619., 605., 606., ..., 1800., 1792., 1727.]], [[2378., 2232., 2686., ..., 2216., 2304., 2038.], [2256., 2289., 2557., ..., 1777., 1998., 2095.], [1933., 1994., 2804., ..., 2027., 2107., 2100.], ..., ... ..., [ 242., 248., 251., ..., 1730., 1616., 1544.], [ 257., 253., 246., ..., 1717., 1608., 1617.], [ 298., 268., 260., ..., 1808., 1774., 1737.]], [[2376., 2035., 2663., ..., 2317., 2828., 2260.], [2161., 2081., 2404., ..., 1831., 2170., 2388.], [1681., 1857., 2790., ..., 2517., 2351., 2390.], ..., [ 201., 199., 212., ..., 1882., 1761., 1696.], [ 225., 202., 186., ..., 1902., 1776., 1788.], [ 242., 228., 214., ..., 2153., 1987., 1938.]], [[2267., 1884., 2269., ..., 2582., 2459., 2093.], [2020., 2059., 2270., ..., 1786., 2236., 2111.], [1536., 1990., 2281., ..., 2023., 2217., 2245.], ..., [ 468., 480., 470., ..., 1726., 1645., 1624.], [ 469., 474., 463., ..., 1694., 1613., 1647.], [ 480., 479., 483., ..., 1803., 1732., 1747.]]], dtype=float32)
- green(time, y, x)float321502.0 1348.0 ... 1328.0 1336.0
- units :
- 1
- crs :
- epsg:32628
- grid_mapping :
- spatial_ref
array([[[1502., 1348., 1825., ..., 1418., 1626., 1270.], [1282., 1418., 1653., ..., 1059., 1355., 1296.], [1019., 1321., 1882., ..., 1391., 1352., 1289.], ..., [ 451., 456., 451., ..., 1201., 1090., 1068.], [ 443., 474., 433., ..., 1092., 1039., 1078.], [ 454., 474., 465., ..., 1088., 1095., 1074.]], [[1627., 1491., 1865., ..., 1493., 1678., 1488.], [1395., 1559., 1750., ..., 1338., 1522., 1544.], [1281., 1601., 1905., ..., 1606., 1540., 1498.], ..., [ 658., 663., 643., ..., 1303., 1228., 1213.], [ 651., 691., 666., ..., 1241., 1219., 1227.], [ 662., 670., 678., ..., 1280., 1248., 1200.]], [[1631., 1534., 1855., ..., 1530., 1667., 1426.], [1528., 1558., 1795., ..., 1280., 1420., 1460.], [1317., 1433., 1944., ..., 1444., 1495., 1444.], ..., ... ..., [ 418., 425., 410., ..., 1321., 1211., 1119.], [ 411., 447., 413., ..., 1308., 1202., 1175.], [ 430., 448., 435., ..., 1376., 1322., 1289.]], [[1610., 1375., 1863., ..., 1632., 2071., 1634.], [1518., 1379., 1647., ..., 1276., 1562., 1701.], [1149., 1306., 1909., ..., 1860., 1683., 1725.], ..., [ 338., 341., 326., ..., 1443., 1301., 1204.], [ 334., 363., 315., ..., 1428., 1297., 1266.], [ 374., 368., 371., ..., 1609., 1487., 1416.]], [[1589., 1391., 1663., ..., 1883., 1867., 1515.], [1479., 1459., 1649., ..., 1345., 1665., 1591.], [1201., 1442., 1599., ..., 1514., 1635., 1658.], ..., [ 616., 633., 621., ..., 1335., 1260., 1220.], [ 616., 625., 595., ..., 1313., 1256., 1241.], [ 622., 634., 633., ..., 1408., 1328., 1336.]]], dtype=float32)
- blue(time, y, x)float32815.0 742.0 ... 1016.0 1002.0
- units :
- 1
- crs :
- epsg:32628
- grid_mapping :
- spatial_ref
array([[[ 815., 742., 1010., ..., 926., 1092., 853.], [ 693., 788., 917., ..., 691., 862., 889.], [ 558., 746., 1035., ..., 903., 883., 864.], ..., [ 230., 255., 248., ..., 838., 732., 734.], [ 240., 251., 224., ..., 748., 707., 741.], [ 240., 239., 251., ..., 739., 727., 712.]], [[ 953., 871., 1098., ..., 1002., 1098., 1006.], [ 818., 919., 1044., ..., 898., 1052., 1028.], [ 769., 984., 1129., ..., 1119., 1039., 1030.], ..., [ 403., 414., 413., ..., 911., 839., 830.], [ 419., 418., 410., ..., 867., 827., 839.], [ 429., 426., 410., ..., 881., 862., 831.]], [[ 946., 915., 1136., ..., 1043., 1131., 974.], [ 921., 925., 1084., ..., 875., 982., 997.], [ 801., 883., 1162., ..., 976., 1025., 1016.], ..., ... ..., [ 191., 186., 191., ..., 989., 889., 792.], [ 191., 195., 167., ..., 990., 864., 839.], [ 211., 182., 186., ..., 1020., 965., 918.]], [[ 777., 654., 970., ..., 1031., 1375., 1121.], [ 674., 693., 804., ..., 803., 995., 1163.], [ 520., 676., 1013., ..., 1250., 1138., 1177.], ..., [ 74., 66., 58., ..., 1022., 893., 800.], [ 68., 61., 47., ..., 1047., 888., 871.], [ 84., 73., 65., ..., 1179., 1051., 960.]], [[ 991., 893., 1066., ..., 1366., 1316., 1140.], [ 915., 951., 1047., ..., 975., 1212., 1181.], [ 788., 955., 1028., ..., 1119., 1213., 1198.], ..., [ 478., 484., 490., ..., 1024., 968., 945.], [ 469., 469., 480., ..., 1038., 954., 950.], [ 477., 481., 486., ..., 1080., 1016., 1002.]]], dtype=float32)
- nir(time, y, x)float323789.0 3519.0 ... 2483.0 2496.0
- units :
- 1
- crs :
- epsg:32628
- grid_mapping :
- spatial_ref
array([[[3789., 3519., 3973., ..., 3178., 3593., 3222.], [3441., 3645., 3866., ..., 2863., 3200., 3165.], [3108., 3637., 4081., ..., 3354., 3288., 3118.], ..., [1871., 1878., 1865., ..., 2795., 2638., 2584.], [1797., 2003., 1787., ..., 2636., 2559., 2647.], [1632., 1877., 1954., ..., 2750., 2729., 2663.]], [[3639., 3498., 3769., ..., 3101., 3402., 3273.], [3346., 3575., 3748., ..., 3110., 3331., 3257.], [3224., 3685., 3798., ..., 3382., 3226., 3224.], ..., [1995., 2030., 1926., ..., 2794., 2718., 2689.], [1891., 2087., 1998., ..., 2747., 2712., 2726.], [1766., 1927., 1972., ..., 2859., 2838., 2739.]], [[3472., 3283., 3723., ..., 3065., 3311., 3157.], [3335., 3405., 3604., ..., 2806., 3068., 3093.], [3079., 3316., 3771., ..., 3093., 3172., 3091.], ..., ... ..., [1962., 2014., 1931., ..., 2377., 2289., 2192.], [1817., 2118., 2058., ..., 2397., 2284., 2274.], [1698., 1895., 2076., ..., 2554., 2459., 2386.]], [[3633., 3515., 4003., ..., 3425., 3979., 3724.], [3450., 3472., 3848., ..., 3340., 3676., 3618.], [3152., 3598., 4054., ..., 3850., 3613., 3646.], ..., [1854., 1923., 1841., ..., 2697., 2558., 2489.], [1769., 2003., 1884., ..., 2708., 2633., 2602.], [1701., 1871., 2007., ..., 3072., 2866., 2748.]], [[3193., 3092., 3455., ..., 3505., 3592., 3149.], [3039., 3184., 3317., ..., 2844., 3264., 3152.], [2710., 3189., 3447., ..., 3312., 3240., 3223.], ..., [1582., 1618., 1626., ..., 2431., 2356., 2356.], [1626., 1685., 1578., ..., 2423., 2374., 2380.], [1546., 1727., 1727., ..., 2574., 2483., 2496.]]], dtype=float32)
- crs :
- epsg:32628
- grid_mapping :
- spatial_ref
Or select a single month:
[5]:
ds.sel(time='2018-05')
[5]:
<xarray.Dataset> Dimensions: (time: 1, x: 1361, y: 1392) Coordinates: * time (time) datetime64[ns] 2018-05-23T11:39:00 * y (y) float64 1.556e+06 1.556e+06 ... 1.528e+06 1.528e+06 * x (x) float64 3.2e+05 3.2e+05 3.201e+05 ... 3.472e+05 3.472e+05 spatial_ref int32 32628 Data variables: red (time, y, x) float32 2626.0 2490.0 2867.0 ... 2050.0 1988.0 green (time, y, x) float32 1895.0 1793.0 2050.0 ... 1584.0 1542.0 blue (time, y, x) float32 1112.0 1080.0 1239.0 ... 1163.0 1132.0 nir (time, y, x) float32 3536.0 3472.0 3682.0 ... 2744.0 2686.0 Attributes: crs: epsg:32628 grid_mapping: spatial_ref
- time: 1
- x: 1361
- y: 1392
- time(time)datetime64[ns]2018-05-23T11:39:00
- units :
- seconds since 1970-01-01 00:00:00
array(['2018-05-23T11:39:00.000000000'], dtype='datetime64[ns]')
- y(y)float641.556e+06 1.556e+06 ... 1.528e+06
- units :
- metre
- resolution :
- -20.0
- crs :
- epsg:32628
array([1555545., 1555525., 1555505., ..., 1527765., 1527745., 1527725.])
- x(x)float643.2e+05 3.2e+05 ... 3.472e+05
- units :
- metre
- resolution :
- 20.0
- crs :
- epsg:32628
array([320025., 320045., 320065., ..., 347185., 347205., 347225.])
- spatial_ref()int3232628
- spatial_ref :
- PROJCS["WGS 84 / UTM zone 28N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-15],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32628"]]
- grid_mapping_name :
- transverse_mercator
array(32628, dtype=int32)
- red(time, y, x)float322626.0 2490.0 ... 2050.0 1988.0
- units :
- 1
- crs :
- epsg:32628
- grid_mapping :
- spatial_ref
array([[[2626., 2490., 2867., ..., 2304., 2370., 2250.], [2573., 2539., 2766., ..., 2149., 2198., 2321.], [2405., 2404., 2935., ..., 2277., 2306., 2264.], ..., [ 677., 698., 693., ..., 2063., 1980., 1917.], [ 679., 675., 683., ..., 2028., 1975., 1970.], [ 730., 702., 697., ..., 2057., 2050., 1988.]]], dtype=float32)
- green(time, y, x)float321895.0 1793.0 ... 1584.0 1542.0
- units :
- 1
- crs :
- epsg:32628
- grid_mapping :
- spatial_ref
array([[[1895., 1793., 2050., ..., 1643., 1719., 1619.], [1854., 1816., 2001., ..., 1540., 1583., 1643.], [1758., 1743., 2119., ..., 1645., 1655., 1629.], ..., [ 718., 737., 709., ..., 1608., 1516., 1482.], [ 684., 705., 709., ..., 1562., 1524., 1502.], [ 702., 707., 694., ..., 1586., 1584., 1542.]]], dtype=float32)
- blue(time, y, x)float321112.0 1080.0 ... 1163.0 1132.0
- units :
- 1
- crs :
- epsg:32628
- grid_mapping :
- spatial_ref
array([[[1112., 1080., 1239., ..., 1133., 1170., 1110.], [1065., 1069., 1192., ..., 1052., 1086., 1147.], [1040., 1065., 1295., ..., 1132., 1157., 1121.], ..., [ 445., 445., 436., ..., 1217., 1125., 1103.], [ 443., 442., 426., ..., 1182., 1134., 1110.], [ 456., 426., 437., ..., 1191., 1163., 1132.]]], dtype=float32)
- nir(time, y, x)float323536.0 3472.0 ... 2744.0 2686.0
- units :
- 1
- crs :
- epsg:32628
- grid_mapping :
- spatial_ref
array([[[3536., 3472., 3682., ..., 3212., 3304., 3284.], [3492., 3494., 3692., ..., 3171., 3232., 3291.], [3357., 3379., 3766., ..., 3206., 3247., 3269.], ..., [1757., 1778., 1679., ..., 2766., 2670., 2614.], [1585., 1771., 1743., ..., 2743., 2690., 2647.], [1520., 1613., 1688., ..., 2765., 2744., 2686.]]], dtype=float32)
- crs :
- epsg:32628
- grid_mapping :
- spatial_ref
Or select a range of dates using slice()
. This selects all observations between the two dates, inclusive of both the start and stop values:
[6]:
ds.sel(time=slice('2018-06', '2019-01'))
[6]:
<xarray.Dataset> Dimensions: (time: 6, x: 1361, y: 1392) Coordinates: * time (time) datetime64[ns] 2018-10-15T11:35:46 ... 2019-01-18T11:... * y (y) float64 1.556e+06 1.556e+06 ... 1.528e+06 1.528e+06 * x (x) float64 3.2e+05 3.2e+05 3.201e+05 ... 3.472e+05 3.472e+05 spatial_ref int32 32628 Data variables: red (time, y, x) float32 1099.0 1051.0 1739.0 ... 1994.0 1951.0 green (time, y, x) float32 1048.0 1058.0 1583.0 ... 1469.0 1409.0 blue (time, y, x) float32 487.0 529.0 836.0 ... 1215.0 1016.0 976.0 nir (time, y, x) float32 3683.0 3957.0 4215.0 ... 2814.0 2792.0 Attributes: crs: epsg:32628 grid_mapping: spatial_ref
- time: 6
- x: 1361
- y: 1392
- time(time)datetime64[ns]2018-10-15T11:35:46 ... 2019-01-...
- units :
- seconds since 1970-01-01 00:00:00
array(['2018-10-15T11:35:46.000000000', '2018-11-19T11:47:30.000000000', '2018-12-09T11:47:28.000000000', '2018-12-14T11:47:31.000000000', '2019-01-03T11:47:35.000000000', '2019-01-18T11:47:33.000000000'], dtype='datetime64[ns]')
- y(y)float641.556e+06 1.556e+06 ... 1.528e+06
- units :
- metre
- resolution :
- -20.0
- crs :
- epsg:32628
array([1555545., 1555525., 1555505., ..., 1527765., 1527745., 1527725.])
- x(x)float643.2e+05 3.2e+05 ... 3.472e+05
- units :
- metre
- resolution :
- 20.0
- crs :
- epsg:32628
array([320025., 320045., 320065., ..., 347185., 347205., 347225.])
- spatial_ref()int3232628
- spatial_ref :
- PROJCS["WGS 84 / UTM zone 28N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-15],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32628"]]
- grid_mapping_name :
- transverse_mercator
array(32628, dtype=int32)
- red(time, y, x)float321099.0 1051.0 ... 1994.0 1951.0
- units :
- 1
- crs :
- epsg:32628
- grid_mapping :
- spatial_ref
array([[[1099., 1051., 1739., ..., 2974., 2897., 2186.], [1075., 1204., 1638., ..., 1696., 2296., 2209.], [ 861., 1314., 1647., ..., 2068., 2182., 2338.], ..., [ 325., 329., 325., ..., 793., 633., 589.], [ 321., 315., 310., ..., 711., 639., 649.], [ 337., 309., 331., ..., 767., 735., 657.]], [[2291., 1949., 2429., ..., 2054., 2651., 2200.], [2159., 1943., 2170., ..., 1784., 1971., 2208.], [1751., 1748., 2549., ..., 2460., 2177., 2185.], ..., [ 242., 248., 251., ..., 1730., 1616., 1544.], [ 257., 253., 246., ..., 1717., 1608., 1617.], [ 298., 268., 260., ..., 1808., 1774., 1737.]], [[2376., 2035., 2663., ..., 2317., 2828., 2260.], [2161., 2081., 2404., ..., 1831., 2170., 2388.], [1681., 1857., 2790., ..., 2517., 2351., 2390.], ..., ... ..., [ 468., 480., 470., ..., 1726., 1645., 1624.], [ 469., 474., 463., ..., 1694., 1613., 1647.], [ 480., 479., 483., ..., 1803., 1732., 1747.]], [[2701., 2189., 2706., ..., 3048., 2902., 2352.], [2520., 2351., 2821., ..., 2030., 2541., 2381.], [1809., 2197., 2761., ..., 2140., 2504., 2546.], ..., [ 333., 336., 350., ..., 1917., 1804., 1786.], [ 329., 330., 333., ..., 1901., 1802., 1808.], [ 362., 339., 339., ..., 2083., 1935., 1947.]], [[2739., 2301., 2918., ..., 2210., 2466., 2147.], [1970., 2465., 2707., ..., 1807., 2402., 2349.], [1763., 2706., 2899., ..., 2415., 2265., 2265.], ..., [ 318., 337., 327., ..., 1861., 1761., 1747.], [ 348., 329., 338., ..., 1917., 1839., 1849.], [ 390., 366., 376., ..., 2203., 1994., 1951.]]], dtype=float32)
- green(time, y, x)float321048.0 1058.0 ... 1469.0 1409.0
- units :
- 1
- crs :
- epsg:32628
- grid_mapping :
- spatial_ref
array([[[1048., 1058., 1583., ..., 2416., 2391., 1738.], [ 985., 1193., 1493., ..., 1453., 1874., 1851.], [ 847., 1231., 1530., ..., 1735., 1827., 1899.], ..., [ 474., 496., 489., ..., 954., 879., 851.], [ 502., 492., 439., ..., 910., 870., 873.], [ 488., 495., 497., ..., 943., 890., 870.]], [[1535., 1352., 1722., ..., 1499., 2018., 1615.], [1519., 1373., 1547., ..., 1323., 1497., 1691.], [1250., 1272., 1759., ..., 1856., 1621., 1664.], ..., [ 418., 425., 410., ..., 1321., 1211., 1119.], [ 411., 447., 413., ..., 1308., 1202., 1175.], [ 430., 448., 435., ..., 1376., 1322., 1289.]], [[1610., 1375., 1863., ..., 1632., 2071., 1634.], [1518., 1379., 1647., ..., 1276., 1562., 1701.], [1149., 1306., 1909., ..., 1860., 1683., 1725.], ..., ... ..., [ 616., 633., 621., ..., 1335., 1260., 1220.], [ 616., 625., 595., ..., 1313., 1256., 1241.], [ 622., 634., 633., ..., 1408., 1328., 1336.]], [[1879., 1533., 1937., ..., 2213., 2127., 1683.], [1709., 1636., 1997., ..., 1401., 1840., 1698.], [1280., 1534., 1900., ..., 1574., 1814., 1814.], ..., [ 470., 473., 477., ..., 1448., 1343., 1299.], [ 476., 491., 453., ..., 1430., 1324., 1322.], [ 486., 490., 476., ..., 1566., 1448., 1429.]], [[1893., 1567., 2051., ..., 1541., 1741., 1497.], [1415., 1714., 1881., ..., 1295., 1723., 1651.], [1277., 1882., 2022., ..., 1713., 1572., 1573.], ..., [ 443., 463., 436., ..., 1379., 1288., 1233.], [ 455., 470., 458., ..., 1427., 1346., 1322.], [ 486., 481., 485., ..., 1655., 1469., 1409.]]], dtype=float32)
- blue(time, y, x)float32487.0 529.0 836.0 ... 1016.0 976.0
- units :
- 1
- crs :
- epsg:32628
- grid_mapping :
- spatial_ref
array([[[ 487., 529., 836., ..., 1800., 1777., 1288.], [ 461., 566., 793., ..., 993., 1356., 1393.], [ 405., 639., 783., ..., 1222., 1318., 1455.], ..., [ 235., 221., 212., ..., 518., 421., 400.], [ 196., 214., 211., ..., 462., 440., 424.], [ 231., 215., 224., ..., 497., 482., 438.]], [[ 815., 686., 968., ..., 1022., 1429., 1165.], [ 767., 722., 835., ..., 884., 1037., 1196.], [ 645., 703., 982., ..., 1309., 1139., 1192.], ..., [ 191., 186., 191., ..., 989., 889., 792.], [ 191., 195., 167., ..., 990., 864., 839.], [ 211., 182., 186., ..., 1020., 965., 918.]], [[ 777., 654., 970., ..., 1031., 1375., 1121.], [ 674., 693., 804., ..., 803., 995., 1163.], [ 520., 676., 1013., ..., 1250., 1138., 1177.], ..., ... ..., [ 478., 484., 490., ..., 1024., 968., 945.], [ 469., 469., 480., ..., 1038., 954., 950.], [ 477., 481., 486., ..., 1080., 1016., 1002.]], [[1023., 854., 1106., ..., 1518., 1454., 1180.], [ 916., 894., 1146., ..., 964., 1271., 1213.], [ 686., 876., 1080., ..., 1057., 1258., 1279.], ..., [ 241., 265., 265., ..., 1099., 989., 934.], [ 231., 254., 250., ..., 1076., 979., 952.], [ 266., 266., 248., ..., 1159., 1054., 1058.]], [[ 941., 849., 1136., ..., 1023., 1158., 990.], [ 704., 912., 1028., ..., 833., 1139., 1103.], [ 683., 1060., 1124., ..., 1142., 1052., 1068.], ..., [ 197., 184., 194., ..., 983., 905., 853.], [ 201., 193., 193., ..., 1033., 934., 916.], [ 226., 220., 214., ..., 1215., 1016., 976.]]], dtype=float32)
- nir(time, y, x)float323683.0 3957.0 ... 2814.0 2792.0
- units :
- 1
- crs :
- epsg:32628
- grid_mapping :
- spatial_ref
array([[[3683., 3957., 4215., ..., 4262., 4238., 3548.], [3490., 3989., 4297., ..., 3533., 3687., 3582.], [3317., 3860., 4371., ..., 3788., 3658., 3712.], ..., [2067., 2045., 2055., ..., 2970., 3016., 3089.], [2203., 2205., 1973., ..., 3020., 3036., 2979.], [1907., 2230., 2175., ..., 3105., 3026., 3096.]], [[3124., 3258., 3661., ..., 3070., 3690., 3385.], [3125., 3076., 3589., ..., 3139., 3432., 3317.], [2917., 3316., 3666., ..., 3579., 3313., 3312.], ..., [1962., 2014., 1931., ..., 2377., 2289., 2192.], [1817., 2118., 2058., ..., 2397., 2284., 2274.], [1698., 1895., 2076., ..., 2554., 2459., 2386.]], [[3633., 3515., 4003., ..., 3425., 3979., 3724.], [3450., 3472., 3848., ..., 3340., 3676., 3618.], [3152., 3598., 4054., ..., 3850., 3613., 3646.], ..., ... ..., [1582., 1618., 1626., ..., 2431., 2356., 2356.], [1626., 1685., 1578., ..., 2423., 2374., 2380.], [1546., 1727., 1727., ..., 2574., 2483., 2496.]], [[3867., 3648., 4274., ..., 4106., 4300., 3688.], [3773., 3751., 4053., ..., 3208., 3784., 3660.], [3205., 3678., 4147., ..., 3764., 3791., 3748.], ..., [1893., 1911., 1915., ..., 2724., 2597., 2627.], [1961., 2075., 1857., ..., 2725., 2635., 2638.], [1860., 2040., 2064., ..., 2908., 2794., 2809.]], [[3840., 3777., 4013., ..., 3285., 3707., 3486.], [3330., 3838., 4009., ..., 3341., 3643., 3548.], [3425., 3929., 4012., ..., 3567., 3486., 3491.], ..., [1949., 1946., 1794., ..., 2678., 2572., 2561.], [1866., 2004., 1924., ..., 2746., 2677., 2661.], [1792., 1961., 1900., ..., 3056., 2814., 2792.]]], dtype=float32)
- crs :
- epsg:32628
- grid_mapping :
- spatial_ref
To select the nearest time to a desired time value, we set it to use a nearest neighbour method, 'nearest'
. We have to specify the time using a datetime
object, otherwise xarray indexing assumes we are selecting a range, like the ds.sel(time='2018-05')
month example above.
Here, we have picked a date at the start of December 2018. 'nearest'
will find the observation closest to that date.
[7]:
target_time = np.datetime64('2018-12-01')
ds.sel(time=target_time, method='nearest')
[7]:
<xarray.Dataset> Dimensions: (x: 1361, y: 1392) Coordinates: time datetime64[ns] 2018-12-09T11:47:28 * y (y) float64 1.556e+06 1.556e+06 ... 1.528e+06 1.528e+06 * x (x) float64 3.2e+05 3.2e+05 3.201e+05 ... 3.472e+05 3.472e+05 spatial_ref int32 32628 Data variables: red (y, x) float32 2376.0 2035.0 2663.0 ... 2153.0 1987.0 1938.0 green (y, x) float32 1610.0 1375.0 1863.0 ... 1609.0 1487.0 1416.0 blue (y, x) float32 777.0 654.0 970.0 867.0 ... 1179.0 1051.0 960.0 nir (y, x) float32 3633.0 3515.0 4003.0 ... 3072.0 2866.0 2748.0 Attributes: crs: epsg:32628 grid_mapping: spatial_ref
- x: 1361
- y: 1392
- time()datetime64[ns]2018-12-09T11:47:28
- units :
- seconds since 1970-01-01 00:00:00
array('2018-12-09T11:47:28.000000000', dtype='datetime64[ns]')
- y(y)float641.556e+06 1.556e+06 ... 1.528e+06
- units :
- metre
- resolution :
- -20.0
- crs :
- epsg:32628
array([1555545., 1555525., 1555505., ..., 1527765., 1527745., 1527725.])
- x(x)float643.2e+05 3.2e+05 ... 3.472e+05
- units :
- metre
- resolution :
- 20.0
- crs :
- epsg:32628
array([320025., 320045., 320065., ..., 347185., 347205., 347225.])
- spatial_ref()int3232628
- spatial_ref :
- PROJCS["WGS 84 / UTM zone 28N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-15],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32628"]]
- grid_mapping_name :
- transverse_mercator
array(32628, dtype=int32)
- red(y, x)float322376.0 2035.0 ... 1987.0 1938.0
- units :
- 1
- crs :
- epsg:32628
- grid_mapping :
- spatial_ref
array([[2376., 2035., 2663., ..., 2317., 2828., 2260.], [2161., 2081., 2404., ..., 1831., 2170., 2388.], [1681., 1857., 2790., ..., 2517., 2351., 2390.], ..., [ 201., 199., 212., ..., 1882., 1761., 1696.], [ 225., 202., 186., ..., 1902., 1776., 1788.], [ 242., 228., 214., ..., 2153., 1987., 1938.]], dtype=float32)
- green(y, x)float321610.0 1375.0 ... 1487.0 1416.0
- units :
- 1
- crs :
- epsg:32628
- grid_mapping :
- spatial_ref
array([[1610., 1375., 1863., ..., 1632., 2071., 1634.], [1518., 1379., 1647., ..., 1276., 1562., 1701.], [1149., 1306., 1909., ..., 1860., 1683., 1725.], ..., [ 338., 341., 326., ..., 1443., 1301., 1204.], [ 334., 363., 315., ..., 1428., 1297., 1266.], [ 374., 368., 371., ..., 1609., 1487., 1416.]], dtype=float32)
- blue(y, x)float32777.0 654.0 970.0 ... 1051.0 960.0
- units :
- 1
- crs :
- epsg:32628
- grid_mapping :
- spatial_ref
array([[ 777., 654., 970., ..., 1031., 1375., 1121.], [ 674., 693., 804., ..., 803., 995., 1163.], [ 520., 676., 1013., ..., 1250., 1138., 1177.], ..., [ 74., 66., 58., ..., 1022., 893., 800.], [ 68., 61., 47., ..., 1047., 888., 871.], [ 84., 73., 65., ..., 1179., 1051., 960.]], dtype=float32)
- nir(y, x)float323633.0 3515.0 ... 2866.0 2748.0
- units :
- 1
- crs :
- epsg:32628
- grid_mapping :
- spatial_ref
array([[3633., 3515., 4003., ..., 3425., 3979., 3724.], [3450., 3472., 3848., ..., 3340., 3676., 3618.], [3152., 3598., 4054., ..., 3850., 3613., 3646.], ..., [1854., 1923., 1841., ..., 2697., 2558., 2489.], [1769., 2003., 1884., ..., 2708., 2633., 2602.], [1701., 1871., 2007., ..., 3072., 2866., 2748.]], dtype=float32)
- crs :
- epsg:32628
- grid_mapping :
- spatial_ref
You can select the closest time before a given time using ffill
(forward-fill).
[8]:
previous_time = ds.sel(time=target_time, method='ffill')
previous_time.blue.plot();

To select the closest time after a given time, use bfill
(back-fill).
[9]:
next_time = ds.sel(time=target_time, method='bfill')
next_time.blue.plot()
[9]:
<matplotlib.collections.QuadMesh at 0x7f07d004ff60>

The same methods also work on a list of times:
[10]:
many_times = np.array([
'2018-06-23',
'2018-09-13',
'2018-11-02'
], dtype=np.datetime64)
nearest = ds.sel(time=many_times, method='nearest')
nearest.blue.plot(col='time', vmin=0);

Using the datetime accessor¶
xarray
allows you to easily extract additional information from the time
dimension in Digital Earth Africa data. For example, we can get a list of what season each observation belongs to:
[11]:
ds.time.dt.season
[11]:
<xarray.DataArray 'season' (time: 34)> array(['DJF', 'DJF', 'DJF', 'DJF', 'MAM', 'MAM', 'MAM', 'MAM', 'SON', 'SON', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'JJA', 'JJA', 'JJA', 'SON', 'SON', 'SON', 'SON', 'SON', 'DJF'], dtype='<U3') Coordinates: * time (time) datetime64[ns] 2018-01-08T11:46:56 ... 2019-12-19T11:... spatial_ref int32 32628
- time: 34
- 'DJF' 'DJF' 'DJF' 'DJF' 'MAM' 'MAM' ... 'SON' 'SON' 'SON' 'SON' 'DJF'
array(['DJF', 'DJF', 'DJF', 'DJF', 'MAM', 'MAM', 'MAM', 'MAM', 'SON', 'SON', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'DJF', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'MAM', 'JJA', 'JJA', 'JJA', 'SON', 'SON', 'SON', 'SON', 'SON', 'DJF'], dtype='<U3')
- time(time)datetime64[ns]2018-01-08T11:46:56 ... 2019-12-...
- units :
- seconds since 1970-01-01 00:00:00
array(['2018-01-08T11:46:56.000000000', '2018-01-23T11:41:48.000000000', '2018-02-17T11:37:13.000000000', '2018-02-22T11:42:49.000000000', '2018-03-04T11:47:19.000000000', '2018-03-09T11:46:52.000000000', '2018-04-03T11:47:15.000000000', '2018-05-23T11:39:00.000000000', '2018-10-15T11:35:46.000000000', '2018-11-19T11:47:30.000000000', '2018-12-09T11:47:28.000000000', '2018-12-14T11:47:31.000000000', '2019-01-03T11:47:35.000000000', '2019-01-18T11:47:33.000000000', '2019-02-12T11:47:37.000000000', '2019-02-22T11:55:41.000000000', '2019-03-09T11:47:38.000000000', '2019-03-14T12:00:04.000000000', '2019-03-24T11:59:54.000000000', '2019-04-08T11:59:55.000000000', '2019-04-18T11:48:41.000000000', '2019-05-03T11:47:45.000000000', '2019-05-13T11:47:45.000000000', '2019-05-23T11:47:46.000000000', '2019-05-28T11:47:40.000000000', '2019-06-07T11:47:39.000000000', '2019-07-02T11:47:45.000000000', '2019-07-07T11:47:42.000000000', '2019-09-30T11:47:37.000000000', '2019-10-20T11:47:38.000000000', '2019-10-30T11:47:38.000000000', '2019-11-24T11:47:37.000000000', '2019-11-29T11:47:32.000000000', '2019-12-19T11:47:32.000000000'], dtype='datetime64[ns]')
- spatial_ref()int3232628
- spatial_ref :
- PROJCS["WGS 84 / UTM zone 28N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-15],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32628"]]
- grid_mapping_name :
- transverse_mercator
array(32628, dtype=int32)
Or the day of the year:
[12]:
ds.time.dt.dayofyear
[12]:
<xarray.DataArray 'dayofyear' (time: 34)> array([ 8, 23, 48, 53, 63, 68, 93, 143, 288, 323, 343, 348, 3, 18, 43, 53, 68, 73, 83, 98, 108, 123, 133, 143, 148, 158, 183, 188, 273, 293, 303, 328, 333, 353]) Coordinates: * time (time) datetime64[ns] 2018-01-08T11:46:56 ... 2019-12-19T11:... spatial_ref int32 32628
- time: 34
- 8 23 48 53 63 68 93 143 288 ... 158 183 188 273 293 303 328 333 353
array([ 8, 23, 48, 53, 63, 68, 93, 143, 288, 323, 343, 348, 3, 18, 43, 53, 68, 73, 83, 98, 108, 123, 133, 143, 148, 158, 183, 188, 273, 293, 303, 328, 333, 353])
- time(time)datetime64[ns]2018-01-08T11:46:56 ... 2019-12-...
- units :
- seconds since 1970-01-01 00:00:00
array(['2018-01-08T11:46:56.000000000', '2018-01-23T11:41:48.000000000', '2018-02-17T11:37:13.000000000', '2018-02-22T11:42:49.000000000', '2018-03-04T11:47:19.000000000', '2018-03-09T11:46:52.000000000', '2018-04-03T11:47:15.000000000', '2018-05-23T11:39:00.000000000', '2018-10-15T11:35:46.000000000', '2018-11-19T11:47:30.000000000', '2018-12-09T11:47:28.000000000', '2018-12-14T11:47:31.000000000', '2019-01-03T11:47:35.000000000', '2019-01-18T11:47:33.000000000', '2019-02-12T11:47:37.000000000', '2019-02-22T11:55:41.000000000', '2019-03-09T11:47:38.000000000', '2019-03-14T12:00:04.000000000', '2019-03-24T11:59:54.000000000', '2019-04-08T11:59:55.000000000', '2019-04-18T11:48:41.000000000', '2019-05-03T11:47:45.000000000', '2019-05-13T11:47:45.000000000', '2019-05-23T11:47:46.000000000', '2019-05-28T11:47:40.000000000', '2019-06-07T11:47:39.000000000', '2019-07-02T11:47:45.000000000', '2019-07-07T11:47:42.000000000', '2019-09-30T11:47:37.000000000', '2019-10-20T11:47:38.000000000', '2019-10-30T11:47:38.000000000', '2019-11-24T11:47:37.000000000', '2019-11-29T11:47:32.000000000', '2019-12-19T11:47:32.000000000'], dtype='datetime64[ns]')
- spatial_ref()int3232628
- spatial_ref :
- PROJCS["WGS 84 / UTM zone 28N",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-15],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH],AUTHORITY["EPSG","32628"]]
- grid_mapping_name :
- transverse_mercator
array(32628, dtype=int32)
Grouping and resampling by time¶
xarray
also provides some shortcuts for aggregating data over time. In the example below, we first group our data by season, then take the median of each group. This produces a new dataset with only four observations (one per season).
[13]:
# Group the time series into seasons, and take median of each time period
ds_seasonal = ds.groupby('time.season').median(dim='time')
# Plot the output
ds_seasonal.nir.plot(col='season', col_wrap=4)
plt.show()

We can also use the .resample()
method to summarise our dataset into larger chunks of time. In the example below, we produce a median composite for every 6 months of data in our dataset:
[14]:
# Resample to combine each 6 months of data into a median composite
ds_resampled = ds.resample(time="6m").median()
# Plot the new resampled data
ds_resampled.nir.plot(col="time")
plt.show()

Interpolating new timesteps¶
Sometimes, we want to return data for specific times/dates that weren’t observed by a satellite. To estimate what the landscape appeared like on certain dates, we can use the .interp()
method to interpolate between the nearest two observations.
By default, the interp()
method uses linear interpolation (method='linear'
). Another useful option is method='nearest'
, which will return the nearest satellite observation to the specified date(s).
[15]:
# New dates to interpolate data for
new_dates = ['2018-07-25', '2018-09-01', '2018-12-05']
# Interpolate Landsat values for three new dates
ds_interp = ds.interp(time=new_dates)
# Plot the new interpolated data
ds_interp.nir.plot(col='time')
plt.show()

Additional information¶
License: The code in this notebook is licensed under the Apache License, Version 2.0. Digital Earth Africa 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 Slack channel 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.
Compatible datacube version:
[16]:
print(datacube.__version__)
1.8.4.dev63+g6ee0462c
Last Tested:
[17]:
from datetime import datetime
datetime.today().strftime('%Y-%m-%d')
[17]:
'2021-05-20'