Espinhal Cleaned Parcels Analysis

Introduction

This is paper endeavours to shed some light on the use of Normalised Difference Vegetation Index (NDVI) coupled with satellite imagery to check whether areas designated as a Defensible Space) by Portuguese government, has been complied with by clearing the area from trees. One such area, that was presumably cleared in July 2018, is defined by the following polygon north of Espinhal:

Espinhal%20DS.png

In the following analysis we're going to analyse satellite images acquired by Sentinel 2, which will help in calculating the NDVI of this area, and ascertain whether the area under this polygon was indeed cleaned in July 2018.

Initialisation

Import packages required to undergo the analysis

In [34]:
import gdal
import os
import json
import rasterio.plot
import numpy as np
import pandas as pd
import geopandas as gpd
import rasterio as rio
import rasterio.mask as mask
import matplotlib.pyplot as plt
import pandas as pd
from rasterio.plot import show
from rasterio.plot import show_hist
from scipy.stats import ttest_rel

import warnings
warnings.filterwarnings('ignore')

Set working directory

In [ ]:
os.chdir('../Github/DataProjects/Espinhal NDVI Analysis')

Create functions that will be used later. The first function jp2tiff_Wgs84 will convert JP2 images to Tiff, and convert the Coordinate Reference System (CRS) to WGS84, otherwise designated by EPSG number 4326. More on this in later sections.

In [5]:
def jp2tiff_wgs84 (path, dsn):
    #Use GDAL to open JP2 image
    jp2=gdal.Open(path)
    
    #Convert JP2 to Tiff AND convert CRS
    gdal.Warp(dsn, jp2, dstSRS='EPSG:4326')

The following function readandcrop will:

  • Read Tiff
  • Crop tiff based on polygon

More on this in later sections

In [6]:
def readandcrop (path, poly):
    #Open connection with Tiff file
    with rio.open(path) as src:
        
        #Read cropped version of tiff based on Polygon
        raster_image, raster_transform = mask.mask(dataset=src,shapes=[poly],crop=True)
    
    #Return cropped image and its transform
    return raster_image, raster_transform

Shape file: Read, Convert CRS and Generate JSON

Read and Check CRS

As per image above, the shape file which contains the area/polygon that was presumably cleaned on 27/07/2018 will be read:

In [8]:
#Read Shape file
cleaned_parcels=gpd.read_file('shape files/ESF_Flopen_2018.shp')

#Inspect its CRS system:
cleaned_parcels.crs
Out[8]:
{'init': 'epsg:3763'}

It can be seen from above, that the CRS of the shape file is EPSG: 3763, which is also ETRS. As such we're going to convert CRS to EPSG: 4326, also known as WGS84, as this is more widely used and supported:

Convert CRS to WGS84

In [9]:
#Convert to WGS84 (EPSG: 4326)
cleaned_wgs84=cleaned_parcels.to_crs({'init':'epsg:4326'})

#Inspect its CRS system:
cleaned_wgs84.crs
Out[9]:
{'init': 'epsg:4326'}

Extract Polygon and Convert to JSON

Since the shape file read has multiple features/polygons, we're going to extract the polygon based on the cleaning date (27/07/2018), after it'll be converted to JSON, as it is needed for cropping later. Finally, we'll plot the polygon to visually inspect we have the right polygon.

In [10]:
#Extract the parcel in question
cleaned_one_poly=cleaned_wgs84[cleaned_wgs84['DAT_FIM']=='27 de jul*']

#Convert to JSON, which will be used later
cleaned_json=json.loads(cleaned_one_poly.to_json())['features'][0]['geometry']

#Plot for visual validation
cleaned_one_poly.plot()
Out[10]:
<matplotlib.axes._subplots.AxesSubplot at 0x1241beda0>

Read and Convert Rasters

As per above, in this section we're going to read and convert JP2 Sentinel 2 images to Tiffs, as Tiffs are readily more supported and easier to work with. We'll also convert the CRS of the tiff to EPSG: 4326 so as to match the shape file's CRS. In terms of the required file to calculate the NDVI, we are specifically after Band8 (Near Infrared: NIR) and Ban4 (Infrared: IR). For each band, we're going to acquire an image for June and an image for August, specifically the following dates:

  • 2018/06/19
  • 2018/08/18
In [11]:
jp2tiff_wgs84("2018-06/T29TNE_20180619T112111_B04.jp2","B04_1806.tiff")
jp2tiff_wgs84("2018-06/T29TNE_20180619T112111_B08.jp2","B08_1806.tiff")
jp2tiff_wgs84("2018-08/T29TNE_20180818T112111_B04.jp2","B04_1808.tiff")
jp2tiff_wgs84("2018-08/T29TNE_20180818T112111_B08.jp2","B08_1808.tiff")

Read Tiff and Crop to Polygon's Extent

Following on from the previous section, we're going to read a "cropped" version of each tiff produced above, this is to reduce the size of data being used, and makes its manipulation easier and faster. The crop will be defined by the polygon's extent and this will define the crop window used.

In [12]:
b4_06_raster, b4_06_transform=readandcrop("B04_1806.tiff",cleaned_json)
b8_06_raster, b8_06_transform=readandcrop("B08_1806.tiff",cleaned_json)
b4_08_raster, b4_08_transform=readandcrop("B04_1808.tiff",cleaned_json)
b8_08_raster, b8_08_transform=readandcrop("B08_1808.tiff",cleaned_json)

Calculate NDVI

In this section we're going to calculate the NDVI for June and August of the defined area, which will give us an idea of the state of vegetation before and after the assumed cleaning date. If we witness a decrease in NDVI, then there's a strong case that the cleaning has happened. The NDVI formula is defined as follows:

$$NDVI=\frac{NIR - IR}{NIR + IR}$$
In [13]:
ndvi_06=(b8_06_raster - b4_06_raster)/(b8_06_raster + b4_06_raster)
ndvi_08=(b8_08_raster - b4_08_raster)/(b8_08_raster + b4_08_raster)

Plot NDVI for June and August 2018

First we'll inspect visually if there's discernible difference in NDVI before and after the presumed "cleaning" of the area, as defined by the polygon.

In [14]:
#Setup canvas
fig, (ax1, ax2)=plt.subplots(1,2, figsize=(15,15))

#Plot NDVI for June
show(ndvi_06,
     ax=ax1,
     cmap='RdYlGn', 
     vmin=-1,
     vmax=1,
     extent=rio.plot.plotting_extent(b4_06_raster, b4_06_transform),
     transform=b4_06_transform)

ax1.set_title('June-2018', fontsize=16);

#Plot NDVI for August
show(ndvi_08[0],
     ax=ax2,
     cmap='RdYlGn', 
     vmin=-1,
     vmax=1,
     extent=rio.plot.plotting_extent(b4_06_raster, b4_06_transform),
     transform=b4_06_transform)
ax2.set_title('August-2018', fontsize=16);

#Increase space between plots
plt.subplots_adjust(wspace = 0.3)

It can be seen from the above that August had a considerable decrease in NDVI as a whole. But a more formal approach is taken below:

NDVI Histogram Plot (With Mean)

We'll create a frequency plot for the NDVIs in the area covered with mean superimposed to provide is with better view of the distribution of NDVI:

In [33]:
fig, (ax1, ax2)=plt.subplots(2,1, figsize=(10,10))

#Plot histo for June
ax1.hist(ndvi_06.ravel(), bins=np.arange(0.0, 0.8, 0.01));
ax1.set_title('NDVI June 2018 with mean {}'.format(np.nanmean(ndvi_06)), fontsize=16);
ax1.axvline(x=np.nanmean(ndvi_06), 
            color='r',
            linestyle='dashed', 
            linewidth=1)

#Plot histo for June
ax2.hist(ndvi_08.ravel(), bins=np.arange(0.0, 0.8, 0.01));
ax2.set_title('NDVI August 2018 with mean {}'.format(np.nanmean(ndvi_08)), fontsize=16);
ax2.axvline(x=np.nanmean(ndvi_08), 
            color='r',
            linestyle='dashed', 
            linewidth=1)
Out[33]:
<matplotlib.lines.Line2D at 0x125c2a6d8>

It can be noted from above, that the mean dropped significantly in August and the distribution of NDVI deviated away from normal distribution, which was present in June (indicating a natural growth of vegetation, regression to the mean). But to conclude this study, we'll inspect NDVI further by creating a hypothesis test using Student's t-test, for this particular test, we're going to apply a paired test since we're dealing with the same parcel of land but over time:

In [37]:
stat, p = ttest_rel(ndvi_06[0].ravel(), ndvi_08[0].ravel(), nan_policy='omit')
"p-vale of {}".format(p)
Out[37]:
'p-vale of 6.567255777746875e-38'

Conclusion

It can be seen from above that P-Value is considerably small, giving us a strong confidence (95%) that the NDVI has considerably decreased over the course of 2 months, which is a significant indicator that cleaning has been done. As such we can try and investigate further the use of remote sensing for inspecting such activities when coupled with NDVI.