Detecting Wildfires and Hotspots Using Hyperspectral Data

📆 Project Period	December, 2023
📍 GitHub	gitlab.esa.int

Introduction

This code reports an example of wildfire classification using PRISMA hyperspectral data. The analysis carried out through this software is fully described in the following paper:

Spiller, D., Ansalone, L., Amici, S., Piscini, A., and Mathieu, P. P.: ANALYSIS AND DETECTION OF WILDFIRES BY USING PRISMA HYPERSPECTRAL IMAGERY, Int. Arch. Photogramm. Remote Sens. Spatial Inf. Sci., XLIII-B3-2021, 215â222, 2021 https://doi.org/10.5194/isprs-archives-XLIII-B3-2021-215-2021. https://www.int-arch-photogramm-remote-sens-spatial-inf-sci.net/XLIII-B3-2021/215/2021/

Note that in the paper there is also a temperature estimation analysis that is not reported in this notebook.

The PRISMA input image required for the analysis is:

PRS_L1_STD_OFFL_20191227000827_20191227000831_0001

related to an acquisition in Australia on 27 Decemeber 2019 and almost centered around the coordinates lat: 31Â° 29' 52", lon: 151Â° 18' 22". This analysis uses the level 1 image, but can be extended to level 2 (B or D).

and the request from the online catalog is shown in the cell below. Please, register to the online portal prisma.asi.it to access the data and download the reference image.

Please contact the author of this notebook in case you need further information:

dario.spiller@esa.int

dario.spiller@asi.it

dario.spiller@uniroma1.it

Libraries

from shapely import affinity
from shapely.ops import cascaded_union
from shapely.geometry import Point, mapping, Polygon
from mpl_toolkits.basemap import Basemap
from itertools import compress
import model
from model import *
from numpy import save
from numpy.random import seed
from sklearn.model_selection import train_test_split
from sklearn.metrics import confusion_matrix,classification_report,accuracy_score
from keras.utils import np_utils
from sklearn.preprocessing import RobustScaler
from mpl_toolkits.axes_grid1 import make_axes_locatable

import numpy as np
import pandas as pd
import matplotlib as mpl
import geopandas as gpd
import spectral.io.envi as envi
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import tensorflow as tf
import h5py
import glob
import fiona
import pickle
import os

Paths of the input files

sep = os.path.sep
input_shape = 'input_data' + sep + 'new_fields_shp.shp'
input_shape_test = 'input_data' + sep + 'new_fields_shp_test.shp'
path_hdr = 'input_data' + sep + 'BuildLayerStack-ReflThor-Quac-Smoothed-ritaglio1-SVM-class200.hdr'
path_enp = 'input_data' + sep + 'BuildLayerStack-ReflThor-Quac-Smoothed-ritaglio1-SVM-class200'
PRISMA_folder = 'PRS - Australia with fire'

Auxilary functions definition

def extract_product_level(folder, level):
    '''This function searches for the PRISMA file with the required processing level'''

    files = glob.glob(folder + sep + '**' + sep + '*.he5', recursive=True)
    search_output = list(map(lambda x:  level in x, files))

    return next(compress(files, search_output))

def create_database_for_NN(input_shape, rangeX, rangeY):
    '''This function creates the dataset for the 1-D CNN model'''

    fields_shp = gpd.read_file(input_shape)

    pixels_in_fields = []

    for item,shp in enumerate(fields_shp.iterrows()):
        poly = Polygon(shp[1].geometry)
        class_id = shp[1].CLASS_ID
        class_name = shp[1].CLASS_NAME
        pixels_in_fields.append({'id':item,'CLASS_ID':class_id,'CLASS_NAME':class_name,'pixels':[]})
        check = 0
        print('\\r', item + 1, '/', len(fields_shp),  end='    ')
        for i in rangeX: #range(lon.shape[0]):
            for j in rangeY: #range(lon.shape[1]):
                p = Point(lon[i,j],lat[i,j])
                if p.within(poly):
                    pixels_in_fields[item]['pixels'].append(DATA_CUBE[i,:,j].T)
                    check = 1
        if check == 0:
            print('Attenzione!')

    return pixels_in_fields

def create_database_for_final_map(rangeX, rangeY):
    '''This function creates the image for the final segmentation map'''

    N_elem = (lon[rangeX.start:rangeX.stop, rangeY.start:rangeY.stop]).size
    pixels_in_map = []
    item = 0
    for i in rangeX: #range(lon.shape[0]):
        for j in rangeY: #range(lon.shape[1]):
            p = Point(lon[i,j],lat[i,j])
            pixels_in_map.append({'id':item,'CLASS_ID':np.random.randint(5),'CLASS_NAME':'unknown','pixels':[DATA_CUBE[i,:,j].T]})
            item += 1
        print('\\r', f'{item} / {N_elem}',  end= '    ')

    return pixels_in_map

def save_data_pickle(input_data, name):
    '''This function save the data by using pickle'''

    data = []
    labels = []

    for element in input_data:
        for row in element['pixels']:
            data.append(np.array(row).reshape(len(row),1))
            labels.append(element['CLASS_ID'])

    data = np.array(data)
    labels = np.array(labels)

    pippo = (data,labels)

    output = open(f'model/data/{name}.pkl', 'wb')
    pickle.dump(pippo, output)
    output.close()

    return 'Done!'

Ground truth data

Training data

fields_shp = gpd.read_file(input_shape)
fields_shp.iloc[:13]

Test data

fields_shp = gpd.read_file(input_shape_test)
fields_shp.iloc[:13]

print('The classes considered for the classification are:')
for class_ in np.unique(fields_shp.CLASS_NAME):
    print(' - ' + class_)

The classes considered for the classification are:
 - Baresoil
 - Burned
 - Fire
 - Smoke
 - Vegetation

Classification map from SVM analysis

The ground truth reference data are related to the analysis reported in the following paper:

Amici, S.; Piscini, A. Exploring PRISMA Scene for Fire Detection: Case Study of 2019 Bushfires in Ben Halls Gap National Park, NSW, Australia. Remote Sens. 2021, 13, 1410. https://doi.org/10.3390/rs13081410

img = envi.open(path_hdr, path_enp)
print(img)

print('Metadata associated to the reference SVM classification:\\n')
metadata = img.metadata
for key, value in metadata.items():
    print(' - ', key, ' : ', value)

plt.imshow(np.squeeze(img.asarray()))
plt.gcf().set_size_inches(18, 9)
plt.show()

PRISMA Reference image

print('Loading PRISMA image...')
Level_1  = extract_product_level(PRISMA_folder,'L1')
ref_path = os.getcwd() + sep + Level_1
f = h5py.File(ref_path,'r')
f.visit(print)

print(f.attrs.keys())

# SWIR Channel SWIR: 920 â 2505 nm
SWIR = f['HDFEOS']['SWATHS'][f'PRS_{Level_1.split("_")[1]}_HCO']['Data Fields']['SWIR_Cube']
print(SWIR)
# VNIR Channel VNIR: 400 â 1010 nm
VNIR = f['HDFEOS']['SWATHS'][f'PRS_{Level_1.split("_")[1]}_HCO']['Data Fields']['VNIR_Cube']
print(VNIR)

VNIR_array = VNIR[()][:,3:,:] # deleting the all-zero elements
SWIR_array = SWIR[()][:,:-2,:]

def search_band(bands,cw,target_wavelength):
    if len(cw) != bands.shape[1]:
        raise ValueError('Dimensions do not match!')
    else:
        temp = np.abs(cw - target_wavelength)
        target_band_index = np.where(temp == temp.min())[0]
        print(f'Index of target wavelength {target_wavelength} nm is {target_band_index}')
        target_band = np.squeeze(bands[:,target_band_index,:])
        return target_band

SWIR_CW = f.attrs.get('List_Cw_Swir')[:-2]
VNIR_CW = f.attrs.get('List_Cw_Vnir')[3:]

if len(VNIR_CW) != VNIR_array.shape[1]:
    raise ValueError('Dimensions do not match!')
else:
    IMG_660 = search_band(VNIR_array,VNIR_CW,660) #660
    IMG_760 = search_band(VNIR_array,VNIR_CW,760)
    IMG_770 = search_band(VNIR_array,VNIR_CW,770)
    IMG_780 = search_band(VNIR_array,VNIR_CW,780)


if len(SWIR_CW) != SWIR_array.shape[1]:
    raise ValueError('Dimensions do not match!')
else:
    IMG_1100 = search_band(SWIR_array,SWIR_CW,1100) #1100
    IMG_1700 = search_band(SWIR_array,SWIR_CW,1700) #1700
    IMG_1990 = search_band(SWIR_array,SWIR_CW,1990)
    IMG_2010 = search_band(SWIR_array,SWIR_CW,2010)
    IMG_2040 = search_band(SWIR_array,SWIR_CW,2040)
    IMG_2060 = search_band(SWIR_array,SWIR_CW,2060)
    IMG_2430 = search_band(SWIR_array,SWIR_CW,2430)
    IMG_1088 = search_band(SWIR_array,SWIR_CW,1088)
    IMG_2200 = search_band(SWIR_array,SWIR_CW,2200)

lat = f[f'HDFEOS/SWATHS/PRS_{Level_1.split("_")[1]}_HCO/Geolocation Fields/Latitude_VNIR'][:]
lon = f[f'HDFEOS/SWATHS/PRS_{Level_1.split("_")[1]}_HCO/Geolocation Fields/Longitude_VNIR'][:]

import matplotlib.pyplot as plt
import numpy as np

red = VNIR_array[:,35,:] # 36=632.13165 nm, 35=641.33325 nm
green = VNIR_array[:,45,:] # 48=530.66705 nm, 45=554.5646 nm
blue = VNIR_array[:,57,:] # 463.731 nm
IM = np.array([red,green,blue])#.reshape(1181,1203,3)
w,h = IM.shape[1],IM.shape[2]
print('Dimensions of the image:', w,h)
IM = IM.reshape(3,w,h)
print('Dimensions of the pre-processed RGB image:',IM.shape)

t=(w,h,3)
A=np.zeros(t,dtype=np.uint16)
for i in range(w):
    for j in range(h):
        A[i,j]=[red[i,j],green[i,j],blue[i,j]]
print('Dimensions of the post-processed RGB image:',A.shape)

# Array containing both SWIR and VNIR
DATA_CUBE = np.concatenate((SWIR_array, VNIR_array), axis=1)

RGB Composite

# initialize the figure
fig,ax = plt.subplots(figsize = (16,16))

# Define lats and longs for the figure
lons_fig = [np.min(lon), np.max(lon)]
lats_fig = [np.min(lat), np.max(lat)]
lat_0 = np.mean(lats_fig)
lon_0 = np.mean(lons_fig)

# initialize the map
map_kwargs = dict(projection='merc', resolution='l',
                  llcrnrlat=np.min(lats_fig), urcrnrlat=np.max(lats_fig),
                  llcrnrlon=np.min(lons_fig), urcrnrlon=np.max(lons_fig),
                  lat_0=lat_0, lon_0=lon_0)
m = Basemap(**map_kwargs)

# drawing parallels and meridians
m.drawparallels(np.arange(np.min(lats_fig),np.max(lats_fig),0.2*(np.max(lats_fig)-np.min(lats_fig))),labels=[1,0,0,0],
                color=[0.6,0.6,0.6], fontsize=30, rotation=30, fmt='%8.5g')
m.drawmeridians(np.arange(np.min(lons_fig),np.max(lons_fig),0.2*(np.max(lons_fig)-np.min(lons_fig))),labels=[0,0,0,1],
                color=[0.6,0.6,0.6], fontsize=30, rotation=30, fmt='%8.5g')

# drawing the polygons in the shapefile
m.readshapefile(os.getcwd() + sep + input_shape[:-4], 'shapefile', drawbounds = False)

for info,shape in zip(m.shapefile_info,m.shapefile):
    if 'Baresoil' in info['CLASS_NAME']:
        line_color = 'pink'
    elif 'Burned' in info['CLASS_NAME']:
        line_color = 'brown'
    elif 'Fire' in info['CLASS_NAME']:
        line_color = 'y'
    elif 'Smoke' in info['CLASS_NAME']:
        line_color = 'blue'
    elif 'Vegetation' in info['CLASS_NAME']:
        line_color = 'green'
    x, y = zip(*shape)
    m.plot(x, y, marker=None, color=line_color, linewidth=1)

# drawing the PRISMA image
mesh_rgb = A[:, :-1, :]/np.max(A)  # also A[:, 1:, :] works fine
colorTuple = mesh_rgb.reshape((mesh_rgb.shape[0] * mesh_rgb.shape[1]), 3)
alpha = 1.0
colorTuple = np.insert(colorTuple,3,alpha,axis=1)

## What you put in for the image doesn't matter because of the color mapping
m.pcolormesh(lon, lat, A[:,:,0], latlon=True, color=colorTuple)

# Showing image
plt.show()

PRISMA images with overlapped ground-truth points

# initialize the figure
fig,ax = plt.subplots(figsize = (16,16))

# Define lats and longs for the figure
lons_fig = [151.2, 151.26]
lats_fig = [-31.629,-31.590]
lat_0 = np.mean(lats_fig)
lon_0 = np.mean(lons_fig)

# initialize the map
map_kwargs = dict(projection='merc', resolution='l',
                  llcrnrlat=np.min(lats_fig), urcrnrlat=np.max(lats_fig),
                  llcrnrlon=np.min(lons_fig), urcrnrlon=np.max(lons_fig),
                  lat_0=lat_0, lon_0=lon_0)
m = Basemap(**map_kwargs)

# drawing parallels and meridians
m.drawparallels(np.arange(np.min(lats_fig),np.max(lats_fig),0.2*(np.max(lats_fig)-np.min(lats_fig))),labels=[1,0,0,0],
                color=[0.6,0.6,0.6], fontsize=30, rotation=30, fmt='%8.5g')
m.drawmeridians(np.arange(np.min(lons_fig),np.max(lons_fig),0.2*(np.max(lons_fig)-np.min(lons_fig))),labels=[0,0,0,1],
                color=[0.6,0.6,0.6], fontsize=30, rotation=30, fmt='%8.5g')

## What you put in for the image doesn't matter because of the color mapping
m.pcolormesh(lon, lat, SWIR_array[:,0,:]/(2**16-1), latlon=True)

# drawing the polygons in the shapefile
m.readshapefile(os.getcwd() + sep + input_shape_test[:-4], 'shapefile1', drawbounds = False)
for info,shape in zip(m.shapefile1_info,m.shapefile1):
    if 'Baresoil' in info['CLASS_NAME']:
        line_color = 'pink'
    elif 'Burned' in info['CLASS_NAME']:
        line_color = 'brown'
    elif 'Fire' in info['CLASS_NAME']:
        line_color = 'y'
    elif 'Smoke' in info['CLASS_NAME']:
        line_color = 'blue'
    elif 'Veg' in info['CLASS_NAME']:
        line_color = 'green'
    x, y = zip(*shape)
    m.plot(x, y, marker=None, color=line_color, linewidth=3)

m.readshapefile(os.getcwd() + sep + input_shape[:-4], 'shapefile2', drawbounds = False)
for info,shape in zip(m.shapefile2_info,m.shapefile2):
    if 'Baresoil' in info['CLASS_NAME']:
        line_color = 'pink'
    elif 'Burned' in info['CLASS_NAME']:
        line_color = 'brown'
    elif 'Fire' in info['CLASS_NAME']:
        line_color = 'y'
    elif 'Smoke' in info['CLASS_NAME']:
        line_color = 'blue'
    elif 'Veg' in info['CLASS_NAME']:
        line_color = 'green'
    x, y = zip(*shape)
    m.plot(x, y, marker=None, color=line_color, linewidth=1)


# create an axes on the right side of ax. The width of cax will be 5%
# of ax and the padding between cax and ax will be fixed at 0.05 inch.
divider = make_axes_locatable(ax)
cax = divider.append_axes("right", size="5%", pad=0.05)

cb = plt.colorbar(cax=cax)
cb.ax.tick_params(labelsize=25)

# Showing image
plt.show()

# initialize the figure
fig,ax = plt.subplots(figsize = (16,16))

# Define lats and longs for the figure
lons_fig = [np.min(lon), np.max(lon)]
lats_fig = [np.min(lat), np.max(lat)]
lat_0 = np.mean(lats_fig)
lon_0 = np.mean(lons_fig)

# initialize the map
map_kwargs = dict(projection='merc', resolution='l',
                  llcrnrlat=np.min(lats_fig), urcrnrlat=np.max(lats_fig),
                  llcrnrlon=np.min(lons_fig), urcrnrlon=np.max(lons_fig),
                  lat_0=lat_0, lon_0=lon_0)
m = Basemap(**map_kwargs)

# drawing parallels and meridians
m.drawparallels(np.arange(np.min(lats_fig),np.max(lats_fig),0.2*(np.max(lats_fig)-np.min(lats_fig))),labels=[1,0,0,0],
                color=[0.6,0.6,0.6], fontsize=30, rotation=30, fmt='%8.5g')
m.drawmeridians(np.arange(np.min(lons_fig),np.max(lons_fig),0.2*(np.max(lons_fig)-np.min(lons_fig))),labels=[0,0,0,1],
                color=[0.6,0.6,0.6], fontsize=30, rotation=30, fmt='%8.5g')

# drawing the polygons in the shapefile
m.readshapefile(os.getcwd() + sep + input_shape[:-4], 'shapefile', drawbounds = False)

for info,shape in zip(m.shapefile_info,m.shapefile):
    if 'Baresoil' in info['CLASS_NAME']:
        line_color = 'pink'
    elif 'Burned' in info['CLASS_NAME']:
        line_color = 'brown'
    elif 'Fire' in info['CLASS_NAME']:
        line_color = 'y'
    elif 'Smoke' in info['CLASS_NAME']:
        line_color = 'blue'
    elif 'Vegetation' in info['CLASS_NAME']:
        line_color = 'green'
    x, y = zip(*shape)
    m.plot(x, y, marker=None, color=line_color, linewidth=1)

# drawing the PRISMA image
mesh_rgb = A[550:700, 850:999, :]/np.max(A)  # also A[:, 1:, :] works fine
colorTuple = mesh_rgb.reshape((mesh_rgb.shape[0] * mesh_rgb.shape[1]), 3)
alpha = 1.0
colorTuple = np.insert(colorTuple,3,alpha,axis=1)

## What you put in for the image doesn't matter because of the color mapping
m.pcolormesh(lon[550:700, 850:1000], lat[550:700, 850:1000], A[550:700,850:1000,0], latlon=True, color=colorTuple)
plt.show()

Extraction of the datasets for the CNN model

# TEST data
rangeX = range(500,800)
rangeY = range(700,1000)

pixels_in_fields_test = create_database_for_NN(input_shape_test, rangeX, rangeY)
save_data_pickle(pixels_in_fields_test,'data_test')

# TRAIN data
rangeX = range(500,800)
rangeY = range(700,1000)

pixels_in_fields = create_database_for_NN(input_shape, rangeX, rangeY)
save_data_pickle(pixels_in_fields,'data')

# SOUTH MAP data
rangeX = range(550,700)
rangeY = range(850,1000)

pixels_in_map = create_database_for_final_map(rangeX, rangeY)
save_data_pickle(pixels_in_map,'data_for_map_south')

# NORTH MAP data
rangeX = range(300,600)
rangeY = range(150,450)

pixels_in_map = create_database_for_final_map(rangeX, rangeY)
save_data_pickle(pixels_in_map,'data_for_map_north')

Classification with 1-D CNN model

seed(1)
tf.random.set_seed(2)

gpus = tf.config.experimental.list_physical_devices('GPU')
if gpus:
    try:
        # Currently, memory growth needs to be the same across GPUs
        for gpu in gpus:
            #tf.config.experimental.set_virtual_device_configuration(gpus[0],
                #[tf.config.experimental.VirtualDeviceConfiguration(memory_limit=1000)])
            tf.config.experimental.set_memory_growth(gpu, True)
        logical_gpus = tf.config.experimental.list_logical_devices('GPU')
        print(len(gpus), "Physical GPUs,", len(logical_gpus), "Logical GPUs")
    except RuntimeError as e:
        # Memory growth must be set before GPUs have been initialized
        print(e)


ref_path = 'model/data'
train_path = 'model/data'
model_path = 'model/model_PRISMA.h5'

# Loading the dataset
pkl_file = open(ref_path + sep + 'data.pkl', 'rb')
data1 = pickle.load(pkl_file)
pkl_file.close()
data0 = data1[0]
data = data0/np.max(data0)
labels = data1[1]

# Loading the dataset for the test
pkl_file = open(ref_path + sep + 'data_test.pkl', 'rb')
data1 = pickle.load(pkl_file)
pkl_file.close()
X_test0 = data1[0]
X_test = X_test0/np.max(X_test0)
y_test0 = data1[1] - 1 # il ' - 1' serve a far partire i labels da zero!

pkl_file = open(ref_path + sep + 'data_for_map_south.pkl', 'rb')
data1 = pickle.load(pkl_file)
pkl_file.close()
X_map_south0 = data1[0]
X_map_south = X_map_south0/np.max(X_test0)
y_test_map_south = data1[1] #

pkl_file = open(ref_path + sep + 'data_for_map_north.pkl', 'rb')
data1 = pickle.load(pkl_file)
pkl_file.close()
X_map_north0 = data1[0]
X_map_north = X_map_north0/np.max(X_test0)
y_test_map_north = data1[1] #


print('Defining training, validation, and test datasets...')
X_train = data
y_train0 = labels - 1 # il ' - 1' serve a far partire i labels da zero!

# Convert labels to categorical variables
y_train = np_utils.to_categorical(y_train0)
y_test = np_utils.to_categorical(y_test0)
y_test_map_south = np_utils.to_categorical(y_test_map_south)
y_test_map_north = np_utils.to_categorical(y_test_map_north)


train_dataset = tf.data.Dataset.from_tensor_slices((X_train,y_train))
test_dataset = tf.data.Dataset.from_tensor_slices((X_test,y_test))

map_south = tf.data.Dataset.from_tensor_slices((X_map_south,y_test_map_south))
map_north = tf.data.Dataset.from_tensor_slices((X_map_north,y_test_map_north))


BATCH_SIZE = 5
SHUFFLE_BUFFER_SIZE = 100

train_batches = train_dataset.shuffle(SHUFFLE_BUFFER_SIZE).batch(BATCH_SIZE)
test_batches = test_dataset.batch(BATCH_SIZE)
map_south = map_south.batch(BATCH_SIZE)
map_north = map_north.batch(BATCH_SIZE)

L = len(X_train[0])
print(L)

# Model construction
model = model_1D(dim=L, learn_rate=1e-4, lmbda=1e-5, drop=0.0, FL=3, n_filters=112,
                  init='he_normal') #  glorot_normal   he_normal

checkpoint = ModelCheckpoint(model_path, save_best_only=True, save_weights_only=True)
stop = EarlyStopping(verbose = 1, patience = 10)
cal = [checkpoint,stop]

print('Fitting the model...')
results = model.fit(train_batches,
                    epochs=200,
                    callbacks=cal,
                    validation_data=test_batches,
                    shuffle='batch')


pred = model.predict(test_batches)
pred_classes = np.argmax(pred, axis=1)
save(ref_path + sep + 'PRISMA_fire_predictions.npy',pred_classes) # SPECIFICARE
print(classification_report(y_test0,pred_classes))
print(accuracy_score(y_test0, pred_classes))
print(confusion_matrix(y_test0, pred_classes))


pred_map_south = model.predict(map_south)
pred_classes = np.argmax(pred_map_south, axis=1)
save(ref_path + sep + 'PRISMA_fire_predictions_map_south.npy',pred_classes) # SPECIFICARE

pred_map_north = model.predict(map_north)
pred_classes = np.argmax(pred_map_north, axis=1)
save(ref_path + sep + 'PRISMA_fire_predictions_map_north.npy',pred_classes) # SPECIFICARE

ref_path = 'model' + sep + 'data' + sep
data = np.load(ref_path + 'PRISMA_fire_predictions_map_south.npy')
data1 = data.reshape(150,150)
colors = ['yellow', 'blue', 'brown', 'green', 'pink']

cmap = mpl.colors.ListedColormap(colors)
plt.imshow(np.fliplr(data1.T)[35:,10:], interpolation='none', cmap=cmap)
plt.axis('off')
plt.savefig("wildfire_south.png", bbox_inches='tight',dpi=300)
plt.gcf().set_size_inches(18, 9)
plt.show()

pippo = img.asarray()

pippo = pippo - 1
print(np.unique(pippo))

for i in range(pippo.shape[0]):
    for j in range(pippo.shape[1]):
        if pippo[i,j]>2:
            pippo[i,j] -= 1

print(np.unique(pippo))

colors = ['pink', 'green', 'brown', 'blue', 'yellow']

cmap = mpl.colors.ListedColormap(colors)
plt.imshow(np.squeeze(pippo), interpolation='none', cmap=cmap)
plt.axis('off')
plt.savefig("wildfire_south_stefania.png", bbox_inches='tight',dpi=300)
plt.gcf().set_size_inches(18, 9)
plt.show()

IDs = {'Fire': 1, 'Smoke': 2, 'Burned': 3, 'Vegetation': 4, 'Baresoil': 5}

ref_path = 'model' + sep + 'data' + sep
data = np.load(ref_path + 'PRISMA_fire_predictions_map_north.npy')
data1 = data.reshape(300,300)
colors = ['yellow', 'blue', 'brown', 'green', 'pink']

cmap = mpl.colors.ListedColormap(colors)
plt.imshow(np.fliplr(data1.T), interpolation='none', cmap=cmap)
plt.axis('off')
# create a patch (proxy artist) for every color
patches = [ mpatches.Patch(facecolor=colors[i],edgecolor='black', label=j) for i,j in enumerate(IDs.keys()) ]
# put those patched as legend-handles into the legend
plt.legend(handles=patches, bbox_to_anchor=(0.19, 1), borderaxespad=0. )


plt.savefig("wildfire_north.png", bbox_inches='tight',dpi=300)
plt.gcf().set_size_inches(18, 9)
plt.show()