Fractional Cover Percentiles

Background: Fractional Cover (FC25) data is stored in DEA. More information on Fractional Cover can be found in the ‘Introduction to Fractional Cover’ notebook.

What does this notebook do?: This notebook demonstrates using Fractional Cover Percentiles (calculated using datacube stats).The data is then saved to NetCDF and GeoTIFF for analysis elsewhere.

Tags: Landsats 5,7,8, products, Fractional Cover, dc.list_products, dc.load, query, plot, image, Plot Multiple Images, quantiles, WOFLs, WOfS, WOfS filtration, NetCDF,:index:GeoTIFF,:index:datacube-stats, Percentile

Import modules from standard libraries, datacube and files

[3]:
%matplotlib inline

import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt
import numpy as np
import os
import sys

#modules for datacube
import datacube
from datacube.storage import masking
from datacube.storage.storage import write_dataset_to_netcdf
from datacube_stats.statistics import Percentile
from datacube.helpers import ga_pq_fuser, write_geotiff
from digitalearthau.utils import wofs_fuser ## this is the group_by for wofs!
# Import external functions from dea-notebooks --note Scripts directory has moved recently
sys.path.append('../10_Scripts/')
import DEAPlotting, DEADataHandling

#ignore datacube warnings (needs to be last import statement)
import warnings
warnings.filterwarnings('ignore', module='datacube')

# set datacube alias (just a string with what you're doing)
dc = datacube.Datacube(app='dc-FC')

Specify the query with the coordinates and date range

[2]:
#the query is a dictionary where the keys are the spatio-temporal specs
query = {
        'lat': (-35.25, -35.35),
        'lon': (149.05, 149.17),
        'time':('2015-01-01', '2017-01-31')
        }

load the data according to our query

[5]:
ds = DEADataHandling.load_clearlandsat(dc, query,product='fc')
Loading ls5 PQ
    Skipping ls5
Loading ls7 PQ
    Loading 0 filtered ls7 timesteps
Loading ls8 PQ
    Loading 15 filtered ls8 timesteps
Combining and sorting ls5, ls7 and ls8 data

Plot Fractional Cover and Unmixing Error Bands

Note that Unmixing Error is high over water

[6]:
#just a scene we want to plot - this can be done other ways
scene = 0
[7]:
#set up our images on a grid using gridspec
plt.figure(figsize=(12,8))
gs = gridspec.GridSpec(2,2) # set up a 2 x 2 grid of 4 images for better presentation

ax1=plt.subplot(gs[0,0])
ds.PV.isel(time=scene).plot(cmap='gist_earth_r')
ax1.set_title('PV')

ax2=plt.subplot(gs[1,0])
ds.BS.isel(time=scene).plot(cmap='Oranges')
ax2.set_title('BS')

ax3=plt.subplot(gs[0,1])
ds.NPV.isel(time=scene).plot(cmap='copper')
ax3.set_title('NPV')

ax4=plt.subplot(gs[1,1])
ds.UE.isel(time=scene).plot(cmap='magma')
ax4.set_title('UE')

plt.tight_layout()
plt.show()
../../_images/notebooks_09_Workflows_Fractional_Cover_Percentiles_10_0.png

load the wofs feature layers (wofls) within the same query as Fractional cover, using ‘like’

[8]:
#load the wofs feature layers (wofls) within the same query as Fractional cover, using 'like'
wofls = dc.load(product = 'wofs_albers', like=ds, group_by='solar_day', fuse_func=wofs_fuser)

Remove areas of water from the Fractional Cover dataset

[28]:
wetwofl = masking.make_mask(wofls, wet=True)
#match WOFL times to our fractional cover times
unwofld = ds.where(ds.time == wetwofl.time)
#mask out water from fractional cover
unwofld = unwofld.where(wetwofl.water==False, drop=False)
del wofls
[11]:
unwofld.isel(time=0).PV.plot()
[11]:
<matplotlib.collections.QuadMesh at 0x7f113e316358>
../../_images/notebooks_09_Workflows_Fractional_Cover_Percentiles_15_1.png

Calculate percentiles of Fractional Cover

This can be done on the command line using datacube-stats, but this notebook demonstrates how to do it in a notebook. The Percentile function can be used on any dataset, but requires you to create an object and then ‘compute’ it.

Specify percentiles to calculate

[30]:
#edit this list if you want different quantiles/percentiles
percentiles = [0, 5, 20, 50, 80, 95, 100]

Calculate percentiles

[32]:
# drop unmixing error and data percentage, which make no sense as quantiles
unwofld =unwofld.drop(['UE', 'data_perc'])
[33]:
unwofld
[33]:
<xarray.Dataset>
Dimensions:  (time: 15, x: 492, y: 500)
Coordinates:
  * time     (time) datetime64[ns] 2015-01-02T23:50:19.500000 ...
  * y        (y) float64 -3.953e+06 -3.953e+06 -3.953e+06 -3.953e+06 ...
  * x        (x) float64 1.542e+06 1.542e+06 1.542e+06 1.542e+06 1.542e+06 ...
Data variables:
    BS       (time, y, x) float64 18.0 25.0 31.0 28.0 32.0 31.0 27.0 26.0 ...
    PV       (time, y, x) float64 45.0 45.0 46.0 53.0 50.0 48.0 54.0 51.0 ...
    NPV      (time, y, x) float64 36.0 29.0 22.0 19.0 17.0 21.0 19.0 21.0 ...
Attributes:
    crs:      EPSG:3577

create Percentile object

[35]:
FC_percents =Percentile(percentiles)

compute percentiles

[36]:
FC_percentiles = FC_percents.compute(unwofld)

reattach attributes

[45]:
#reattatch attributes here
FC_percentiles.attrs = unwofld.attrs
[46]:
FC_percentiles
[46]:
<xarray.Dataset>
Dimensions:     (x: 492, y: 500)
Coordinates:
  * y           (y) float64 -3.953e+06 -3.953e+06 -3.953e+06 -3.953e+06 ...
  * x           (x) float64 1.542e+06 1.542e+06 1.542e+06 1.542e+06 ...
Data variables:
    BS_PC_0     (y, x) float64 14.0 21.0 31.0 28.0 21.0 26.0 23.0 19.0 13.0 ...
    PV_PC_0     (y, x) float64 41.0 37.0 36.0 36.0 31.0 30.0 37.0 39.0 44.0 ...
    NPV_PC_0    (y, x) float64 29.0 19.0 8.0 5.0 8.0 4.0 7.0 15.0 20.0 17.0 ...
    BS_PC_5     (y, x) float64 14.0 24.0 33.0 30.0 27.0 31.0 26.0 24.0 14.0 ...
    PV_PC_5     (y, x) float64 42.0 39.0 39.0 37.0 34.0 33.0 38.0 40.0 45.0 ...
    NPV_PC_5    (y, x) float64 30.0 20.0 8.0 6.0 16.0 4.0 8.0 16.0 23.0 20.0 ...
    BS_PC_20    (y, x) float64 14.0 27.0 38.0 33.0 32.0 36.0 30.0 30.0 16.0 ...
    PV_PC_20    (y, x) float64 42.0 41.0 40.0 40.0 38.0 34.0 43.0 44.0 46.0 ...
    NPV_PC_20   (y, x) float64 33.0 22.0 9.0 12.0 16.0 6.0 9.0 17.0 25.0 ...
    BS_PC_50    (y, x) float64 18.0 30.0 45.0 42.0 36.0 51.0 38.0 32.0 21.0 ...
    PV_PC_50    (y, x) float64 46.0 42.0 42.0 43.0 39.0 40.0 48.0 45.0 48.0 ...
    NPV_PC_50   (y, x) float64 36.0 26.0 12.0 17.0 23.0 12.0 14.0 20.0 29.0 ...
    BS_PC_80    (y, x) float64 20.0 33.0 47.0 43.0 40.0 54.0 46.0 36.0 26.0 ...
    PV_PC_80    (y, x) float64 48.0 45.0 46.0 50.0 48.0 46.0 54.0 51.0 53.0 ...
    NPV_PC_80   (y, x) float64 38.0 29.0 17.0 19.0 27.0 19.0 18.0 23.0 33.0 ...
    BS_PC_95    (y, x) float64 24.0 36.0 50.0 45.0 43.0 62.0 50.0 44.0 29.0 ...
    PV_PC_95    (y, x) float64 50.0 47.0 48.0 51.0 49.0 48.0 55.0 52.0 53.0 ...
    NPV_PC_95   (y, x) float64 38.0 30.0 19.0 22.0 27.0 25.0 20.0 25.0 34.0 ...
    BS_PC_100   (y, x) float64 27.0 43.0 51.0 50.0 53.0 64.0 52.0 44.0 29.0 ...
    PV_PC_100   (y, x) float64 51.0 48.0 49.0 53.0 50.0 49.0 57.0 52.0 54.0 ...
    NPV_PC_100  (y, x) float64 39.0 31.0 22.0 29.0 30.0 27.0 25.0 28.0 39.0 ...
Attributes:
    crs:      EPSG:3577

Plot percentiles

[43]:
FC_percentiles.PV_PC_5.plot(cmap='YlGn', vmax=100)
plt.show()
FC_percentiles.NPV_PC_5.plot(cmap='copper', vmax=100)
plt.show()
FC_percentiles.BS_PC_5.plot(cmap='Oranges_r', vmax=100)
plt.show()
../../_images/notebooks_09_Workflows_Fractional_Cover_Percentiles_30_0.png
../../_images/notebooks_09_Workflows_Fractional_Cover_Percentiles_30_1.png
../../_images/notebooks_09_Workflows_Fractional_Cover_Percentiles_30_2.png

Edit save file path to choose where to save your output files

[54]:
#save files to your home directory (add your outpath here if you want to change it)
savefilepath = os.path.expanduser('~/')

Save Fractional Cover percentiles to NetCDF

  • add a unit to the unitless dimension so that it writes to NetCDF

[51]:
FC_percentiles.attrs['units'] = 'fractional_cover_percentage_percentile'
[52]:
FC_percentiles
[52]:
<xarray.Dataset>
Dimensions:     (x: 492, y: 500)
Coordinates:
  * y           (y) float64 -3.953e+06 -3.953e+06 -3.953e+06 -3.953e+06 ...
  * x           (x) float64 1.542e+06 1.542e+06 1.542e+06 1.542e+06 ...
Data variables:
    BS_PC_0     (y, x) float64 14.0 21.0 31.0 28.0 21.0 26.0 23.0 19.0 13.0 ...
    PV_PC_0     (y, x) float64 41.0 37.0 36.0 36.0 31.0 30.0 37.0 39.0 44.0 ...
    NPV_PC_0    (y, x) float64 29.0 19.0 8.0 5.0 8.0 4.0 7.0 15.0 20.0 17.0 ...
    BS_PC_5     (y, x) float64 14.0 24.0 33.0 30.0 27.0 31.0 26.0 24.0 14.0 ...
    PV_PC_5     (y, x) float64 42.0 39.0 39.0 37.0 34.0 33.0 38.0 40.0 45.0 ...
    NPV_PC_5    (y, x) float64 30.0 20.0 8.0 6.0 16.0 4.0 8.0 16.0 23.0 20.0 ...
    BS_PC_20    (y, x) float64 14.0 27.0 38.0 33.0 32.0 36.0 30.0 30.0 16.0 ...
    PV_PC_20    (y, x) float64 42.0 41.0 40.0 40.0 38.0 34.0 43.0 44.0 46.0 ...
    NPV_PC_20   (y, x) float64 33.0 22.0 9.0 12.0 16.0 6.0 9.0 17.0 25.0 ...
    BS_PC_50    (y, x) float64 18.0 30.0 45.0 42.0 36.0 51.0 38.0 32.0 21.0 ...
    PV_PC_50    (y, x) float64 46.0 42.0 42.0 43.0 39.0 40.0 48.0 45.0 48.0 ...
    NPV_PC_50   (y, x) float64 36.0 26.0 12.0 17.0 23.0 12.0 14.0 20.0 29.0 ...
    BS_PC_80    (y, x) float64 20.0 33.0 47.0 43.0 40.0 54.0 46.0 36.0 26.0 ...
    PV_PC_80    (y, x) float64 48.0 45.0 46.0 50.0 48.0 46.0 54.0 51.0 53.0 ...
    NPV_PC_80   (y, x) float64 38.0 29.0 17.0 19.0 27.0 19.0 18.0 23.0 33.0 ...
    BS_PC_95    (y, x) float64 24.0 36.0 50.0 45.0 43.0 62.0 50.0 44.0 29.0 ...
    PV_PC_95    (y, x) float64 50.0 47.0 48.0 51.0 49.0 48.0 55.0 52.0 53.0 ...
    NPV_PC_95   (y, x) float64 38.0 30.0 19.0 22.0 27.0 25.0 20.0 25.0 34.0 ...
    BS_PC_100   (y, x) float64 27.0 43.0 51.0 50.0 53.0 64.0 52.0 44.0 29.0 ...
    PV_PC_100   (y, x) float64 51.0 48.0 49.0 53.0 50.0 49.0 57.0 52.0 54.0 ...
    NPV_PC_100  (y, x) float64 39.0 31.0 22.0 29.0 30.0 27.0 25.0 28.0 39.0 ...
Attributes:
    crs:      EPSG:3577
    units:    fractional_cover_percentage_percentile

#Write out percentiles to NetCDF

[56]:
try:
    DEADataHandling.write_your_netcdf(FC_percentiles, 'FC_percentiles', savefilepath+'FC_percentiles.nc', crs = ds.crs)
# #complain if the file already exists but don't fail
except RuntimeError as err:
    print("RuntimeError: {0}".format(err))
print('wrote to netCDF' )
wrote to netCDF

Save Fractional Cover percentiles to GeoTiff

[60]:
try:
    ds = FC_percentiles
    write_geotiff(savefilepath+'FC_percentiles.tif', ds)
    print('wrote to GeoTiff' )
#         DEADataHandling.write_your_netcdf(FC_quantiles.isel(quantile=quant), 'FC_Q'+str(quantiles[quant]), savefilepath+'FC_Q_'+(str(quantiles[quant]).replace('.','_'))+'.nc', crs = ds.crs)
except RuntimeError as err:
    print("RuntimeError: {0}".format(err))
wrote to GeoTiff

References

  1. GA, 2015. Fractional Cover (FC25) Product Description. https://d28rz98at9flks.cloudfront.net/79676/Fractional_Cover_FC25_v1_5.PDF

  2. TERN Auscover, 2012. Fractional cover - Landsat, Joint Remote Sensing Research Program algorithm, Australia coverage. Last modified by Peter Scarth on 2017/03/24 09:53. http://data.auscover.org.au/xwiki/bin/view/Product+pages/Landsat+Fractional+Cover