In [2]:
import os
from pyrsgis import raster
from sklearn import cluster
from osgeo import gdal
import geopandas as gdp
import pandas as pd
import numpy as np
import rasterio as rio
import earthpy as et
import earthpy.spatial as es
import earthpy.plot as ep
import matplotlib.pyplot as plt
from glob import glob
from rasterio.plot import show
import pandas as pd
import cv2
Group Members:
Luis Garcia
Caleb Roudy
Rowan Beswick
Introduction¶
**The 2019 to 2020 Australian bushfire season, also known as Black summer, was one of the worst natural disasters in modern history. Not only causing havic within some on the wider ecosystem in central- southern Australia, but also affecting human society. By the end of it all, more than 9,000 buildings were vanquished, along with 34 unfortunate lives lost.These damages can be assessed from above with the many commericial domestic and international satelittes and platforms. One of the more prominent platforms from above is Landsat 8, which has been around monitoring most portions of the globe. The images from the satelitte are utilized to monitor and detect surface processes such as wildfire, ice melt, deforestation, as well as others.
Methods¶
Data used entails (band 1 — band 7) of Landsat 8/OLI Sydney fire as features and try to predict the binary fire/nonfire class. These 2 images will be used for training and testing. Finally, another multispectral Landsat 8 data acquired in the the outside of Newcastle( rural town of Murrurundi) will be used for new predictions. To further prove the effiency o this method, a California fire image will also. Lastly, a Desicion Tree method will also be implemented on the California Fire.
Portions of the electromagnetic spectrum that Landsat 8/OLI covers
This is a graph showing the different ranges and disributions of the electromagnetic spectrum at which each sensors covers. Most areas at which they cover were designed to capture solar radiation signal from the so called "Atmosphereic Window"
Source : https://prd-wret.s3.us-west-2.amazonaws.com/assets/palladium/production/s3fs-public/styles/full_width/public/thumbnails/image/dmidS2LS7Comparison.png
In [2]:
# Path of the file we will be using to idenitfy fire
# Image of the fire in Newcastle( rural town of Murrurundi)
landsat_path = glob("LC08_L1TP_090082_20191231_20200111_01_T1/LC08_L1TP_090082_20191231_20200111_01_T1_B*.TIF")
landsat_path.sort()
array_stack, meta_data = es.stack(landsat_path)
# Image of the fire in the outskirts of Sydney, Australia
landsat_path_2 = glob("LC08_L1TP_090083_20191215_20191226_01_T1/LC08_L1TP_090083_20191215_20191226_01_T1_B*_Refl.TIF")
landsat_path_2.sort()
array_stack_2, meta_data_2 = es.stack(landsat_path_2)
In [3]:
title = ['Aerosol' , 'Blue' , 'Green' , 'Red' , 'NIR' , 'SWIR1' , 'SWIR2']
ep.plot_bands(array_stack,
title=title)
Out[3]:
In [12]:
%matplotlib notebook
# Tr
ep.plot_rgb(array_stack,
rgb=[5, 4, 2],
title="Landsat 8 RGB Image\n of the rural town of Murrurundi (image for neaural network)")
plt.figure(figsize=(4,4))
# Adjust the amount of linear stretch to futher brighten the image
ep.plot_rgb(array_stack_2,
rgb=[5, 4, 2],
title="Landsat 8 RGB Image\n of outside of Sydney (trained image)")
plt.figure(figsize=(4,4))
Up above are images utilizing the portions of the electromagnetic spectrum that Landsat 8 covers. Notice how cloud/smoke cover dissapates as electromagnetic spectrum shifts to longer wavelengths. This process is explained by the scattering proporties of particles in our atmosphere.
The second set of images here are true color composites, using blue, green, red bands of the spectrum
K Means Cluster; finding ideal number of clusters that the data naturally binds to:
In [6]:
# Run the Kmeans algorithm and get the index of data points clusters
file = 'LC08_L1TP_044032_20181108_20181116_01_T1.tar/LC08_L1TP_044032_20181108_20181116_01_T1/LC08_L1TP_044032_20181108_20181116_01_T1_B7.tif'
dataset = gdal.Open(file)
Band= dataset.GetRasterBand(1)
img = Band.ReadAsArray()
X = img.reshape((-1,1))
sse = []
list_k = list(range(1, 5))
for k in list_k:
km =cluster.KMeans(n_clusters=k)
km.fit(X)
sse.append(km.inertia_)
# Plot sse against k
plt.figure(figsize=(6, 6))
plt.plot(list_k, sse, '-o')
plt.xlabel(r'Number of clusters *k*')
plt.ylabel('Sum of squared distance')
Out[6]:
Neural Network Classification¶
This diagram illustrates the theory behind neural netoworks for detecton. In short, features, or bands of an image in this case, are put into the network, as well as a true label of those bands, and ultimately outputs a final map
In [3]:
# Landsat 8/OLI image of Outside of Sydney
southerAus = 'LC08_L1TP_090083_20191215_20191226_01_T1/LC08_20191215_1226_LayersStacked_ALL.tif'
# 'LC08_L1TP_090083_20191215_20191226_01_T1/LC08_20191215_1226_LayersStacked_ALL.tif'
# Binary Classification of fire non fire ground truth of Sydney (produced from K-means)
fire_Binary = 'LC08_L1TP_090083_20191215_20191226_01_T1_NEW/min_distance_fire_nonfire.tif'
# Outside of Newcastle in the rural town of Murrurundi
eastAus = 'LC08_L1TP_090082_20191231_20200111_01_T1/LC08_20191231_20200111_01_T1_Stacked.tif'
#LC08_L1TP_091086_20191222_20200110_01_T1 (Victoria, Australia- outside of the town of Bumberrah)
# Read the rasters as array
ds7, featuresSouthAus = raster.read(southerAus, bands='all')
ds6, labelsAus = raster.read(fire_Binary, bands=1)
ds5, featuresEastAus = raster.read(eastAus, bands='all')
In [13]:
print("Fire b7 image shape: ", featuresSouthAus.shape)
print("Fire b6 image shape: ", labelsAus.shape)
print("Fire b5 image shape: ", featuresEastAus.shape)
In [4]:
from pyrsgis.convert import changeDimension
featuresSouthAus = changeDimension(featuresSouthAus)
labelsAus = changeDimension(labelsAus)
featuresEastAus = changeDimension(featuresEastAus)
nBands = featuresSouthAus.shape[1]
labelsAus= (labelsAus == 1).astype(int)
print("Souther Aus multispectral image shape: ", featuresSouthAus.shape)
print("Southern Aus Fire/nonFire binary shape: ", labelsAus.shape)
print("Victoria multispectral image shape: ", featuresEastAus.shape)
In [6]:
from sklearn.model_selection import train_test_split
#X (Indepent Varibles/Data) #y(Dependent/Labels (Fire/NonFire))
xTrain, xTest, yTrain, yTest = train_test_split(featuresSouthAus, labelsAus, test_size=0.4, random_state=42)
print(xTrain.shape)
print(yTrain.shape)
print(xTest.shape)
print(yTest.shape)
In [7]:
# Reshape the data
xTrain = xTrain.reshape((xTrain.shape[0], 1, xTrain.shape[1]))
xTest = xTest.reshape((xTest.shape[0], 1, xTest.shape[1]))
featuresEastAus = featuresEastAus.reshape((featuresEastAus.shape[0], 1, featuresEastAus.shape[1]))
# Print the shape of reshaped data
print(xTrain.shape, xTest.shape, featuresEastAus.shape)
In [8]:
from tensorflow import keras
# Define the parameters of the model
model = keras.Sequential([
keras.layers.Flatten(input_shape=(1, nBands)),
keras.layers.Dense(14, activation='relu'),
keras.layers.Dense(2, activation='softmax')])
# Define the accuracy metrics and parameters
model.compile(optimizer="adam", loss="sparse_categorical_crossentropy", metrics=["accuracy"])
# Run the model
model.fit(xTrain, yTrain, epochs=2)
Out[8]:
In [15]:
np.unique(xTest)
Out[15]:
In [12]:
yTestPredicted = model.predict(xTest)
Out[12]:
In [13]:
np.unique(yTestPredicted)
Out[13]:
In [20]:
from sklearn.metrics import confusion_matrix, precision_score, recall_score
# Predict for test data
yTestPredicted = model.predict(xTest)
yTestPredicted = yTestPredicted[:,1]
# Calculate and display the error metrics
yTestPredicted = (yTestPredicted>0.5).astype(int)
cMatrix = confusion_matrix(yTest, yTestPredicted)
pScore = precision_score(yTest, yTestPredicted)
rScore = recall_score(yTest, yTestPredicted)
print("Confusion matrix: for 14 nodes\n", cMatrix)
print("\nP-Score: %.3f, R-Score: %.3f" % (pScore, rScore))
In [26]:
landsat_path = glob("LC08_L1TP_090082_20191231_20200111_01_T1/LC08_L1TP_090082_20191231_20200111_01_T1_B*.TIF")
landsat_path.sort()
array_stack, meta_data = es.stack(landsat_path)
print(len(landsat_path))
In [28]:
# Adjust the amount of linear stretch to futher brighten the image
ep.plot_rgb(array_stack,
rgb=[5, 4, 2],
title="Landsat RGB Image\n Linear Stretch Applied")
plt.show()
In [24]:
predicted = model.predict(featuresEastAus)
predicted = predicted[:,1]
#Export raster
prediction = np.reshape(predicted, (ds5.RasterYSize, ds5.RasterXSize))
outFile = 'Eastern Australia Fire Prediction.tif'
raster.export(prediction, ds5, filename=outFile, dtype='float')
Final output¶
In [41]:
# Read File Raster Data as Array using Gdal
img_file = "Eastern Australia Fire Prediction.tif"
gtif = gdal.Open(img_file)
georaster = gtif.ReadAsArray()
# Plot image using matplotlib
plt.figure(figsize=(10,10))
plt.imshow(georaster)
Out[41]:
Demonstration 2¶
In [41]:
# Landsat 8/OLI image of Outside of Sydney
southerAus = 'LC08_L1TP_090083_20191215_20191226_01_T1/LC08_20191215_1226_LayersStacked_ALL.tif'
# 'LC08_L1TP_090083_20191215_20191226_01_T1/LC08_20191215_1226_LayersStacked_ALL.tif'
# Binary Classification of fire non fire ground truth of Sydney (produced from K-means)
fire_Binary = 'LC08_L1TP_090083_20191215_20191226_01_T1_NEW/min_distance_fire_nonfire.tif'
# Outside of Newcastle in the rural town of Murrurundi
cal_fire = 'LC08_L1TP_044032_20181108_20181116_01_T1.tar/LC08_L1TP_044032_20181108_20181116_01_T1/LC08_L1TP_044032_20181108_20181116_01_T1_LayerStacked.tif'
#LC08_L1TP_091086_20191222_20200110_01_T1 (Victoria, Australia- outside of the town of Bumberrah)
# Read the rasters as array
ds1, ftSouthAus = raster.read(southerAus, bands='all')
ds2, labelsAus = raster.read(fire_Binary, bands=1)
ds3, ftCal = raster.read(cal_fire, bands='all')
In [33]:
landsat_path_3 = glob("LC08_L1TP_044032_20181108_20181116_01_T1.tar/LC08_L1TP_044032_20181108_20181116_01_T1/LC08_L1TP_044032_20181108_20181116_01_T1_B*.TIF")
landsat_path_3.sort()
array_stack_3, meta_data_3 = es.stack(landsat_path_3)
ep.plot_rgb(array_stack_3,
rgb=[5, 4, 2],
title="Landsat 8 RGB Image\California Fire(image for neaural network)")
plt.figure(figsize=(4,4))
Out[33]:
In [43]:
print("Fire b7 image shape: ", ftSouthAus.shape)
print("Fire b6 image shape: ", labelsAus.shape)
print("Fire b5 image shape: ", ftCal.shape)
In [44]:
from pyrsgis.convert import changeDimension
ftSouthAus = changeDimension(ftSouthAus)
labelsAus = changeDimension(labelsAus)
ftCal = changeDimension(ftCal)
nBands = ftSouthAus.shape[1]
labelsAus= (labelsAus == 1).astype(int)
print("Souther Aus multispectral image shape: ", ftSouthAus.shape)
print("Southern Aus Fire/nonFire binary shape: ", labelsAus.shape)
print("California FIre multispectral image shape: ", ftCal.shape)
In [45]:
from sklearn.model_selection import train_test_split
#X (Indepent Varibles/Data) #y(Dependent/Labels (Fire/NonFire))
xTrain, xTest, yTrain, yTest = train_test_split(ftSouthAus, labelsAus, test_size=0.4, random_state=42)
print(xTrain.shape)
print(yTrain.shape)
print(xTest.shape)
print(yTest.shape)
In [46]:
# Reshape the data
xTrain = xTrain.reshape((xTrain.shape[0], 1, xTrain.shape[1]))
xTest = xTest.reshape((xTest.shape[0], 1, xTest.shape[1]))
ftCal= ftCal.reshape((ftCal.shape[0], 1, ftCal.shape[1]))
# Print the shape of reshaped data
print(xTrain.shape, xTest.shape, ftCal.shape)
In [47]:
from tensorflow import keras
# Define the parameters of the model
model = keras.Sequential([
keras.layers.Flatten(input_shape=(1, nBands)),
keras.layers.Dense(14, activation='relu'),
keras.layers.Dense(2, activation='softmax')])
# Define the accuracy metrics and parameters
model.compile(optimizer="adam", loss="sparse_categorical_crossentropy", metrics=["accuracy"])
# Run the model
model.fit(xTrain, yTrain, epochs=2)
Out[47]:
In [48]:
from sklearn.metrics import confusion_matrix, precision_score, recall_score
# Predict for test data
yTestPredicted = model.predict(xTest)
yTestPredicted = yTestPredicted[:,1]
# Calculate and display the error metrics
yTestPredicted = (yTestPredicted>0.5).astype(int)
cMatrix = confusion_matrix(yTest, yTestPredicted)
pScore = precision_score(yTest, yTestPredicted)
rScore = recall_score(yTest, yTestPredicted)
print("Confusion matrix: for 14 nodes\n", cMatrix)
print("\nP-Score: %.3f, R-Score: %.3f" % (pScore, rScore))
In [49]:
predicted = model.predict(ftCal)
predicted = predicted[:,1]
#Export raster
prediction = np.reshape(predicted, (ds3.RasterYSize, ds3.RasterXSize))
outFile = 'California Fire Prediction.tif'
raster.export(prediction, ds3, filename=outFile, dtype='float')
Output #2¶
In [59]:
img_file = "California Fire Prediction.tif"
gtif = gdal.Open(img_file)
georaster = gtif.ReadAsArray()
plt.figure(figsize=(7,7))
plt.imshow(georaster)
Out[59]: