Forest Degredation in Vietnam using regression analysis on Landsat Imagery

This notebook:

applies a forest-degredation monitoring method across a time series of landsat imagery
displays a rendering of the computed products,
saves computed products off to disk for further validation.

Motivation

This notebook is inspired by a publication titled Assessment of Forest Degradation in Vietnam Using Landsat Time Series Data authored by Vogelmann et al. You can retrieve a copy from mdpi by following this link.

Algorithmic Profile

This algorithm generates a forest degredation product.
The product is derived from Landsat 7 Collection 1 Tier 2 SR imagery taken from USGS data holdings.
Linear regression is run on an NDVI product, the linear coeffiecient (slope) is used as a proxy for forest degredation (NDVI decrease).

Process

For a select year:

Compute the NDVI across all Landsat acquisitons
Select a time frame of n contiguous acquisitions
run linear regression on time-series stack of each ndvi pixel
capture slope in a lat,lon referenced grid

Flow Diagram

Calculations

NDVI

NDVI is a derived index that correlates well with the existance of vegetation.

$$ NDVI = \frac{(NIR - RED)}{(NIR + RED)}$$



In [1]:

    
def NDVI(dataset):
    return (dataset.nir - dataset.red)/(dataset.nir + dataset.red).rename("NDVI")

Linear Regression

The following code runs regression analysison every pixel.

If it looks messy, that's probably because the underlying regression needs to handle nan values in a very peculiar way.



In [2]:

    
import xarray as xr  
import numpy as np  

def _where_not_nan(arr):
    return np.where(np.isfinite(arr))

def _flatten_shallow(arr):
    return arr.reshape(arr.shape[0] * arr.shape[1])  

def per_pixel_linear_trend(pixel_da: xr.DataArray) -> xr.DataArray:
    time_index_length = len(pixel_da.time)  
    
    ys = _flatten_shallow(pixel_da.values)
    xs = np.array(list(range(time_index_length)))

    not_nan = _where_not_nan(ys)[0].astype(int)

    xs = xs[not_nan]
    ys = ys[not_nan] 

    pf = np.polyfit(xs,ys, 1)
    return xr.DataArray(pf[0])

def trend_product(da: xr.DataArray) -> xr.DataArray:
    stacked = da.stack(allpoints = ['latitude', 'longitude'])
    trend = stacked.groupby('allpoints').apply(per_pixel_linear_trend)
    unstacked = trend.unstack('allpoints')
    return unstacked.rename(dict(allpoints_level_0 = "latitude", allpoints_level_1 = "longitude"))

Cloud Masking



In [3]:

    
import numpy as np

def is_clear_l7(dataset):    
    #Create boolean Masks for clear and water pixels
    clear_pixels = dataset.pixel_qa.values == 2 + 64
    water_pixels = dataset.pixel_qa.values == 4 + 64

    a_clean_mask = np.logical_or(clear_pixels, water_pixels)
    return a_clean_mask

Case study

Spatial Extents



In [4]:

    
# Zanzibar, Tanzania
latitude = (-6.2238, -6.1267)
longitude = (39.2298, 39.2909)
date_range = ('2000-01-01', '2001-12-31')

Display basemap of area



In [5]:

    
from utils.data_cube_utilities.dc_display_map import display_map
display_map(latitude = latitude, longitude = longitude)









    Out[5]:

Load Data

Import datacube



In [6]:

    
import datacube
dc = datacube.Datacube()

Load data



In [7]:

    
data = dc.load( latitude = latitude, 
                longitude = longitude,
                product = "ls7_usgs_sr_scene",
                measurements = ['red', 'nir', 'pixel_qa'])

Create a cloud mask

Unclear pixels will be masked with a nan value. We'll drop pixel_qa from the dataset to preserve memory.



In [8]:

    
mask = is_clear_l7(data)  
data = data.drop(['pixel_qa'])

Calculate NDVI



In [9]:

    
data = NDVI(data)

Filter clouded/occluded NDVI readings



In [10]:

    
data = data.where(mask)
del mask

Run regression



In [11]:

    
from time import time 

t1 = time()
data = trend_product(data)
t2 = time()



In [12]:

    
print(t2 - t1)









    



104.10355734825134

Plot trends below threshold



In [13]:

    
%matplotlib inline 
data.where(data< 0).plot(figsize = (16,11), cmap = 'Reds_r')









    Out[13]:





<matplotlib.collections.QuadMesh at 0x7f4e43c4f9b0>



In [14]:

    
(-data).plot(figsize = (16,11))









    Out[14]:





<matplotlib.collections.QuadMesh at 0x7f4e2fcbeb70>

Cited

Deutscher, Janik & Gutjahr, Karlheinz & Perko, Roland & Raggam, Hannes & Hirschmugl, Manuela & Schardt, Mathias. (2017). Humid Tropical Forest Monitoring with Multi-Temporal L-, C- and X-Band SAR Data. 10.1109/Multi-Temp.2017.8035264.