Logo
  • Home
  • About ESA Φ-lab CIN
  • CIN People
  • Opportunities
  • Projects
  • Φ-talks
  • News
→ THE EUROPEAN SPACE AGENCY
Hyperspectral Data and AI to Classify Crops

Hyperspectral Data and AI to Classify Crops

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

Getting started

To use the code and dataset locally, you'll need a few things:

  1. The code from the repository.
  2. The dataset.
  3. The python environment set up.
  4. Jupyterlab to run the notebook.

Installing the environment

  • conda create -n env_name
  • conda activate env_name
  • conda install ipykernel
  • python -m ipykernel install --user --name=env_name

Importing the required libraries

  • conda install pip
  • pip install -r requirements.txt

Running the notebook locally

  • jupyter lab

Notebook

Processing of PRISMA image using GDAL

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

Processing of PRISMA image using H5PY

print('Loading PRISMA image...')
f = h5py.File(raster_path,'r')
print('List of data in the PRISMA file...')
f.visit(print)
# SWIR Channel SWIR: 920 – 2505 nm
SWIR = f['HDFEOS']['SWATHS']['PRS_L2D_HCO']['Data Fields']['SWIR_Cube']
# VNIR Channel VNIR: 400 – 1010 nm
VNIR = f['HDFEOS']['SWATHS']['PRS_L2D_HCO']['Data Fields']['VNIR_Cube']
print(f.attrs.keys())
Dimensions of the VNIR datacube: (1172, 66, 1299)
Dimensions of the image: 1172 1299
Dimensions of the pre-processed RGB image: (3, 1172, 1299)
Dimensions of the post-processed RGB image: (1172, 1299, 3)

# Importing the shapefile for CASI
shape_path = cwd + os.path.sep + 'acquisitions' + os.path.sep + 'shapefile_EPSG4326.shp'
fields_shp = gpd.read_file(shape_path)
fields_shp
image

Plotting the Tuscany region

png

Plotting the PRISMA image over the Tuscany region

png

RGB composition of the PRISMA image

C:\\Users\\dario\\anaconda3\\lib\\site-packages\\ipykernel_launcher.py:24: MatplotlibDeprecationWarning:
The dedent function was deprecated in Matplotlib 3.1 and will be removed in 3.3. Use inspect.cleandoc instead.
png

Crop fields over the PRISMA image

png

Collecting the PRISMA pixels within the crop fields

Plotting the corn and the tomato spectral curves

png

Creating the training and the test datasets

tomato points:  487
corn points:  1690

Defining the 1D model

Training and testing the model

Logo

About ESA EO

About CIN

About Pi School

ESA Φ-lab Website

ESA Φ-lab Linkedin community

Copyright 2025 @ European Space Agency. All rights reserved.

LinkedInXGitHubInstagramFacebookYouTube
from mpl_toolkits.basemap import Basemap
import matplotlib.patches as mpatches
from matplotlib.patches import Polygon
from matplotlib.collections import PatchCollection

import tensorflow as tf
from tensorflow.keras.layers import Conv1D, MaxPooling1D, Dropout, Flatten, Dense
from tensorflow.keras import Model
from tensorflow.keras.optimizers import Adam
from tensorflow.keras.regularizers import l2
from tensorflow.keras.callbacks import ModelCheckpoint, EarlyStopping
from numpy import save

import os
import os.path

import numpy as np
from numpy.random import seed


from sklearn.model_selection import train_test_split
from sklearn.metrics import confusion_matrix,classification_report

import pickle

import h5py
import geopandas as gpd
import gdal
import matplotlib.pyplot as plt
import shapely
shapely.speedups.disable()

## List input raster files
#raster_path = 'D:\\\\Documenti_Lavoro\\\\PRISMA\\\\data from CNR\\\\PRISMA images\\\\PRS_L2D_STD_20200731101044_20200731101048_0001\\\\PRS_L2D_STD_20200731101044_20200731101048_0001.he5'
cwd = os.getcwd()
raster_path = cwd + os.path.sep + 'PRISMA images' + os.path.sep + 'PRS_L2D_STD_20200731101044_20200731101048_0001' + os.path.sep + 'PRS_L2D_STD_20200731101044_20200731101048_0001.he5'

## Open HDF file
hdflayer = gdal.Open(raster_path, gdal.GA_ReadOnly)

# Get the file metadata
metadata = hdflayer.GetMetadata()
print('Printing all the PRISMA metadata...\\n')
print(metadata)

# List of the datasets within the file
datasets_list = hdflayer.GetSubDatasets()
#print('\\n\\n'.join(map(str, datasets_list)))

# Extract the bands information for the VNIR
SWIR = datasets_list[0][0]
VNIR = datasets_list[2][0]
bands_SWIR = gdal.Open(SWIR, gdal.GA_ReadOnly)
bands_VNIR = gdal.Open(VNIR, gdal.GA_ReadOnly)
VNIR_array = bands_VNIR.ReadAsArray()[:,3:,:] # deleting the all-zero elements
SWIR_array = bands_SWIR.ReadAsArray()[:,:-2,:] # deleting the all-zero elements
band_array = bands_VNIR.ReadAsArray()
<KeysViewHDF5 ['Acquisition_Purpose', 'Acquisition_Size', 'Acquisition_Station', 'Acquisition_Type', 'Atm_LutGeomInfo_RelativeAzimuth', 'Atm_LutGeomInfo_SunZenith', 'Atm_LutGeomInfo_ViewZenith', 'Atm_Lut_version', 'Atmo_RTM_info', 'Atmo_profile_info', 'Aux_SunEarthDistance', 'Aux_SunIrradiance', 'CNM_L2_BINNING', 'CNM_L2_BIN_ON', 'CNM_L2_BSEL_ON', 'CNM_L2_HGRP', 'CNM_PAN_ACQ', 'CNM_SWIR_ACQ', 'CNM_SWIR_SELECT', 'CNM_VNIR_ACQ', 'CNM_VNIR_SELECT', 'Cloudy_pixels_percentage', 'DEM_info', 'Epsg_Code', 'Exit_Code', 'Frame_Type', 'GCP_info', 'ISF_ID_Start', 'Image_ID', 'Integration_Time', 'L1_Processor_Version', 'L1_Quality_CCPerc', 'L2ScalePanMax', 'L2ScalePanMin', 'L2ScaleSwirMax', 'L2ScaleSwirMin', 'L2ScaleVnirMax', 'L2ScaleVnirMin', 'L2d_Quality_flags', 'List_Cw_Swir', 'List_Cw_Swir_Flags', 'List_Cw_Vnir', 'List_Cw_Vnir_Flags', 'List_Fwhm_Swir', 'List_Fwhm_Vnir', 'Main_Electornic_Unit', 'Map_AOT_accuracy', 'Map_WV_accuracy', 'Num_Frames', 'Number_of_ISF', 'PANCorruptedFrameList', 'PAN_ACQ', 'PAN_Corrupted_Frame_Percentage', 'PAN_HGRP', 'PAN_HYP_ACT_RESIDUAL_m', 'PAN_HYP_ALT_RESIDUAL_m', 'PAN_HYP_START_SYNC_FRAME', 'PAN_HYP_START_SYNC_SUBFRAME', 'PAN_HYP_STOP_SYNC_FRAME', 'PAN_HYP_STOP_SYNC_SUBFRAME', 'PAN_Int_Time', 'PE_Gain_SWIR', 'PE_Gain_VNIR', 'Pan_N_Int', 'Pan_Num_Frames', 'Prev_Cdp_File_Name', 'Prev_FKdp_File_Name', 'Prev_Gkdp_File_Name', 'Processing_Level', 'Processing_Station', 'Processing_Time', 'Processor_Name', 'Processor_Version', 'Product_ID', 'Product_LLcorner_easting', 'Product_LLcorner_lat', 'Product_LLcorner_long', 'Product_LLcorner_northing', 'Product_LRcorner_easting', 'Product_LRcorner_lat', 'Product_LRcorner_long', 'Product_LRcorner_northing', 'Product_Name', 'Product_StartTime', 'Product_StopTime', 'Product_ULcorner_easting', 'Product_ULcorner_lat', 'Product_ULcorner_long', 'Product_ULcorner_northing', 'Product_URcorner_easting', 'Product_URcorner_lat', 'Product_URcorner_long', 'Product_URcorner_northing', 'Product_center_easting', 'Product_center_lat', 'Product_center_long', 'Product_center_northing', 'Projection_Id', 'Projection_Name', 'Reference_Ellipsoid', 'SWIRCorruptedFrameList', 'SWIR_BNSTART', 'SWIR_BNSTOP', 'SWIR_Corrupted_Frame_Percentage', 'SWIR_HGRP', 'SWIR_X', 'Sea_pixels_percentage', 'Soi_L0a_EO-EOS', 'Soi_Post_Dark_Calibration_L0aFile', 'Soi_Prev_Dark_Calibration_L0aFile', 'Sun_azimuth_angle', 'Sun_zenith_angle', 'Synch_Time', 'VNIRCorruptedFrameList', 'VNIR_BNSTART', 'VNIR_BNSTOP', 'VNIR_Corrupted_Frame_Percentage', 'VNIR_HGRP', 'VNIR_X']>
VNIR = band_array
VNIR_cube = band_array #np.array(VNIR)
print('Dimensions of the VNIR datacube:',VNIR_cube.shape)
red = VNIR[:,35,:] # 36=632.13165 nm, 35=641.33325 nm
green = VNIR[:,45,:] # 48=530.66705 nm, 45=554.5646 nm
blue = VNIR[:,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)

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

# load lats and longs of the PRISMA image
lat = f['HDFEOS/SWATHS/PRS_L2D_HCO/Geolocation Fields/Latitude'][:]
lon = f['HDFEOS/SWATHS/PRS_L2D_HCO/Geolocation Fields/Longitude'][:]

# Define lats and longs for the figure
lons_fig = [6, 19]
lats_fig = [36.3, 47.3]
lat_0 = np.mean(lats_fig)
lon_0 = np.mean(lons_fig)

# initialize the map
map_kwargs = dict(projection='merc', resolution='h',
                  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.1,0.1,0.1], 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.1,0.1,0.1], fontsize=30, rotation=30, fmt='%8.5g')

m.readshapefile(cwd + os.path.sep + 'Italy_shapefile'+ os.path.sep +'Limiti01012018_g'+ os.path.sep +'Reg01012018_g'+ os.path.sep +'Reg01012018_g_WGS84_EPSG4326',
                'shapefile_EPSG4326', drawbounds = True, linewidth=1)

patches   = []

for info,shape in zip(m.shapefile_EPSG4326_info,m.shapefile_EPSG4326):
    if info['DEN_REG']== 'Toscana':
        patches.append( Polygon(np.array(shape), True) )

ax.add_collection(PatchCollection(patches, facecolor= 'y', edgecolor='k', linewidths=1., zorder=2))

m.drawcoastlines(linewidth=3)
m.drawcountries(linewidth=2)
m.etopo()

# Labels on axes
#plt.xlabel('Longitude (degree)', labelpad=80, fontsize=16)
#plt.ylabel('Latitude (degree)', labelpad=100, fontsize=16)

plt.annotate('ITALY', xy=(0.3333, 0.5), xycoords='axes fraction', fontsize=80, color='black', fontweight='bold')
plt.show()


C:\\Users\\dario\\anaconda3\\lib\\site-packages\\ipykernel_launcher.py:19: MatplotlibDeprecationWarning:
The dedent function was deprecated in Matplotlib 3.1 and will be removed in 3.3. Use inspect.cleandoc instead.
C:\\Users\\dario\\anaconda3\\lib\\site-packages\\ipykernel_launcher.py:28: MatplotlibDeprecationWarning:
The dedent function was deprecated in Matplotlib 3.1 and will be removed in 3.3. Use inspect.cleandoc instead.
C:\\Users\\dario\\anaconda3\\lib\\site-packages\\ipykernel_launcher.py:39: MatplotlibDeprecationWarning:
The dedent function was deprecated in Matplotlib 3.1 and will be removed in 3.3. Use inspect.cleandoc instead.

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

# load lats and longs of the PRISMA image
lat = f['HDFEOS/SWATHS/PRS_L2D_HCO/Geolocation Fields/Latitude'][:]
lon = f['HDFEOS/SWATHS/PRS_L2D_HCO/Geolocation Fields/Longitude'][:]

# Define lats and longs for the figure
lons_fig = [9.6, 12.5]
lats_fig = [42.3, 44.5]
lat_0 = np.mean(lats_fig)
lon_0 = np.mean(lons_fig)

# initialize the map
map_kwargs = dict(projection='merc', resolution='f',
                  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.1,0.1,0.1], 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.1,0.1,0.1], fontsize=30, rotation=30, fmt='%8.5g')

m.readshapefile(cwd + os.path.sep + 'Italy_shapefile'+ os.path.sep +'Limiti01012018_g'+ os.path.sep +'Reg01012018_g'+ os.path.sep +'Reg01012018_g_WGS84_EPSG4326',
                'shapefile_EPSG4326', drawbounds = True, linewidth=1)

patches   = []

for info,shape in zip(m.shapefile_EPSG4326_info,m.shapefile_EPSG4326):
    if info['DEN_REG']== 'Toscana':
        x, y = zip(*shape)
        m.plot(x, y, marker=None,color='y',linewidth=5)

m.drawcoastlines(linewidth=3)
m.drawcountries(linewidth=2)
m.etopo()

# 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
eps_lon = -0.001 # additive offset in degree
eps_lat = -0.0012 # additive offset in degree
m.pcolormesh(lon+eps_lon, lat+eps_lat, A[:,:,0], latlon=True, color=colorTuple)

plt.annotate('TUSCANY', xy=(0.3333, 0.5), xycoords='axes fraction', fontsize=80, color='black', fontweight='bold')

# Labels on axes
#plt.xlabel('Longitude (degree)', labelpad=80, fontsize=16)
#plt.ylabel('Latitude (degree)', labelpad=100, fontsize=16)


# Showing image
plt.show()


C:\\Users\\dario\\anaconda3\\lib\\site-packages\\ipykernel_launcher.py:19: MatplotlibDeprecationWarning:
The dedent function was deprecated in Matplotlib 3.1 and will be removed in 3.3. Use inspect.cleandoc instead.
C:\\Users\\dario\\anaconda3\\lib\\site-packages\\ipykernel_launcher.py:28: MatplotlibDeprecationWarning:
The dedent function was deprecated in Matplotlib 3.1 and will be removed in 3.3. Use inspect.cleandoc instead.
C:\\Users\\dario\\anaconda3\\lib\\site-packages\\ipykernel_launcher.py:38: MatplotlibDeprecationWarning:
The dedent function was deprecated in Matplotlib 3.1 and will be removed in 3.3. Use inspect.cleandoc instead.

from mpl_toolkits.basemap import Basemap
import matplotlib.patches as mpatches

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

# load lats and longs of the PRISMA image
lat = f['HDFEOS/SWATHS/PRS_L2D_HCO/Geolocation Fields/Latitude'][:]
lon = f['HDFEOS/SWATHS/PRS_L2D_HCO/Geolocation Fields/Longitude'][:]

# Define lats and longs for the figure
eps_lon = -0.001 # additive offset in degree
eps_lat = -0.0012 # additive offset in degree
lons_fig = [np.min(lon+eps_lon), np.max(lon+eps_lon)]
lats_fig = [np.min(lat+eps_lat), np.max(lat+eps_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')


m.etopo()

# 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+eps_lon, lat+eps_lat, A[:,:,0], latlon=True, color=colorTuple)

# Labels on axes
#plt.xlabel('Longitude (degree)', labelpad=80, fontsize=16)
#plt.ylabel('Latitude (degree)', labelpad=100, fontsize=16)

# Showing image
plt.show()


from mpl_toolkits.basemap import Basemap
import matplotlib.patches as mpatches

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

# load lats and longs of the PRISMA image
lat = f['HDFEOS/SWATHS/PRS_L2D_HCO/Geolocation Fields/Latitude'][:]
lon = f['HDFEOS/SWATHS/PRS_L2D_HCO/Geolocation Fields/Longitude'][:]

# Define lats and longs for the figure
lons_fig = [11.016, 11.100]
lats_fig = [42.814, 42.886]
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(cwd + os.path.sep + 'acquisitions'+ os.path.sep +'shapefile_EPSG4326','shapefile_EPSG4326',
                drawbounds = False)
for info,shape in zip(m.shapefile_EPSG4326_info,m.shapefile_EPSG4326):
    if info['crop_type']== 'corn':
        line_color = 'm'
    elif info['crop_type']== 'tomato':
        line_color = 'y'
    x, y = zip(*shape)
    m.plot(x, y, marker=None, color=line_color, linewidth=6)

# 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
eps_lon = -0.001 # additive offset in degree
eps_lat = -0.0012 # additive offset in degree
m.pcolormesh(lon+eps_lon, lat+eps_lat, A[:,:,0], latlon=True, color=colorTuple)

# Labels on axes
#plt.xlabel('Longitude (degree)', labelpad=80, fontsize=16)
#plt.ylabel('Latitude (degree)', labelpad=100, fontsize=16)

# polygons legend
corn_fields = mpatches.Patch(color='m', label='Corn fields')
tomato_fields = mpatches.Patch(color='y', label='Tomato fields')
plt.legend(handles=[corn_fields,tomato_fields], fontsize=40)

# Showing image
plt.show()


C:\\Users\\dario\\anaconda3\\lib\\site-packages\\ipykernel_launcher.py:22: MatplotlibDeprecationWarning:
The dedent function was deprecated in Matplotlib 3.1 and will be removed in 3.3. Use inspect.cleandoc instead.
C:\\Users\\dario\\anaconda3\\lib\\site-packages\\ipykernel_launcher.py:32: MatplotlibDeprecationWarning:
The dedent function was deprecated in Matplotlib 3.1 and will be removed in 3.3. Use inspect.cleandoc instead.
C:\\Users\\dario\\anaconda3\\lib\\site-packages\\IPython\\core\\pylabtools.py:132: UserWarning: Creating legend with loc="best" can be slow with large amounts of data.
  fig.canvas.print_figure(bytes_io, **kw)

from shapely.geometry import Point, Polygon

lat = f['HDFEOS/SWATHS/PRS_L2D_HCO/Geolocation Fields/Latitude'][:]
lon = f['HDFEOS/SWATHS/PRS_L2D_HCO/Geolocation Fields/Longitude'][:]

eps_lon = -0.001 # additive offset in degree
eps_lat = -0.0012 # additive offset in degree

lon_corrected = lon + eps_lon
lat_corrected = lat + eps_lat

pixels_in_fields = []

for item,shp in enumerate(fields_shp.iterrows()):
    poly = Polygon(shp[1].geometry)
    crop_type = shp[1].crop_type
    origin = shp[1].origin
    pixels_in_fields.append({'id':item,'crop_type':crop_type,'origin':origin,'pixels':[]})
    print(pixels_in_fields[item])
    for i in range(lon_corrected.shape[0]):
        for j in range(lon_corrected.shape[1]):
            p = Point(lon_corrected[i,j],lat_corrected[i,j])
            if p.within(poly):
                pixels_in_fields[item]['pixels'].append(DATA_CUBE[i,:,j].T)


{'id': 0, 'crop_type': 'tomato', 'origin': 'TOMATO2L', 'pixels': []}
{'id': 1, 'crop_type': 'tomato', 'origin': 'TOMATO3L', 'pixels': []}
{'id': 2, 'crop_type': 'corn', 'origin': 'CORN1', 'pixels': []}
{'id': 3, 'crop_type': 'corn', 'origin': 'CORN2', 'pixels': []}
{'id': 4, 'crop_type': 'corn', 'origin': 'CORN3', 'pixels': []}
{'id': 5, 'crop_type': 'tomato', 'origin': 'TOMATO1L', 'pixels': []}
{'id': 6, 'crop_type': 'corn', 'origin': 'CORN4L', 'pixels': []}

temp = pixels_in_fields

TOMATO_pixels = [temp[0]['pixels'],
                 temp[1]['pixels'],
                 temp[5]['pixels']]

#tomato_array = np.zeros((len(TOMATO_pixels[0]), len(TOMATO_pixels[0][0])))
tomato_array = np.zeros((len(TOMATO_pixels[0]) + len(TOMATO_pixels[1]) + len(TOMATO_pixels[2]), len(TOMATO_pixels[0][0])))
k = 0
for i in range(len(TOMATO_pixels)):
    for elem in TOMATO_pixels[i]:
        tomato_array[k] = elem/(2**16) #np.flip(elem/(2**16))
        k = k + 1
means_tomato = np.mean(tomato_array,axis=0)
stds_tomato = np.std(tomato_array,axis=0)

CORN_pixels = [temp[2]['pixels'],
               temp[3]['pixels'],
               temp[4]['pixels'],
               temp[6]['pixels']]

#corn_array = np.zeros((len(CORN_pixels[0]), len(CORN_pixels[0][0])))
corn_array = np.zeros((len(CORN_pixels[0]) + len(CORN_pixels[1]) + len(CORN_pixels[2]) + len(CORN_pixels[3]), len(CORN_pixels[0][0])))
k = 0
for i in range(len(CORN_pixels)):
    for elem in CORN_pixels[i]:
        corn_array[k] = elem/(2**16) #np.flip(elem/(2**16))
        k = k + 1
means_corn = np.mean(corn_array,axis=0)
stds_corn = np.std(corn_array,axis=0)


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

freq = np.flip(np.append(freq_SWIR,freq_VNIR))
freq = np.append(freq_SWIR,freq_VNIR)

to_del = [171,172,173]
freq1 = np.delete(freq, to_del)
means_tomato1 = np.delete(means_tomato, to_del)
stds_tomato1 = np.delete(stds_tomato, to_del)
means_corn1 = np.delete(means_corn, to_del)
stds_corn1 = np.delete(stds_corn, to_del)

x = np.argsort(freq1)

plt.style.use('ggplot')
plt.rc('font', size=14)          # controls default text sizes
fig, ax = plt.subplots(figsize=(8, 3))
plt.plot(freq1[x],means_tomato1[x], label='Tomato')
plt.fill_between(freq1[x],means_tomato1[x]-3*stds_tomato1[x], means_tomato1[x]+3*stds_tomato1[x], alpha=.3)
plt.plot(freq1[x],means_corn1[x],color='blue', label='Corn')
plt.fill_between(freq1[x],means_corn1[x]-3*stds_corn1[x], means_corn1[x]+3*stds_corn1[x],alpha=.3,color='blue')
plt.ylabel('At-surface\\nReflectance',fontsize=16) #($W m-2 sr-1 nm-1$)
plt.xlabel('Wavelength (nm)',fontsize=16)
ax.set_ylim(0, 0.8)
ax.legend(fontsize=16)
plt.tight_layout()
plt.show()




np.random.seed(10)

data = []
labels = []
data_test = []
labels_test = []

k = 0
for element in pixels_in_fields:
    if element['origin'] == 'TOMATO2L' or element['origin'] == 'CORN3':
        if element['origin'] == 'TOMATO2L':
            label = 0
        if element['origin'] == 'CORN3':
            label = 1
        for row in element['pixels']:
            data_test.append(np.array(row).reshape(len(row),1))
            labels_test.append(label)
    else:
        if element['crop_type'] == 'tomato':
            label = 0
        else:
            label = 1
        for row in element['pixels']:
            data.append(np.array(row).reshape(len(row),1))
            labels.append(label)

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


idx = np.random.permutation(len(labels))
data,labels = data[idx,:,:], labels[idx]
idx = np.random.permutation(len(labels_test))
data_test,labels_test = data_test[idx,:,:], labels_test[idx]

output_data = (data,labels)

output = open('model/data/data.pkl', 'wb')
pickle.dump(output_data, output)
output.close()

output_test_data = (data_test,labels_test)

output = open('model/data/data_test.pkl', 'wb')
pickle.dump(output_test_data, output)
output.close()

tomato_points = 0;
corn_points = 0

for item,shp in enumerate(fields_shp.iterrows()):
    if pixels_in_fields[item]['crop_type'] == 'tomato':
        tomato_points = tomato_points + len(pixels_in_fields[item]['pixels'])
    if pixels_in_fields[item]['crop_type'] == 'corn':
        corn_points = corn_points + len(pixels_in_fields[item]['pixels'])

print('tomato points: ', tomato_points)
print('corn points: ', corn_points)

def model_1D(dim, learn_rate, lmbda, drop, FL, n_filters, init):
    """Parameters:

    dim : image dimension (256)
    learn_rate : learning rate for the optimizer (1e-4)
    lmbda : regularization parameter (1e-5)
    drop : dropout percentage (0.15)
    FL : kernel size (3)
    n_filters: number of filters (112)
    init: kernel initializer (he_normal)

    Output:

    model : unet model"""

    img_input = tf.keras.Input(shape=(dim,1))

    C1 = Conv1D(n_filters,FL,padding='same',activation='relu',kernel_initializer=init,kernel_regularizer=l2(lmbda))(img_input)
    MP1 = MaxPooling1D(pool_size=(2), strides=(2))(C1)

    F = Flatten()(MP1)
    D1 = Dense(units = 128, activation = 'relu')(F)
    D2 = Dense(units = 1, activation = 'sigmoid')(D1)

    model = Model(inputs=img_input,outputs=D2)

    opt = Adam(learning_rate=learn_rate)
    model.compile(loss='mean_squared_error',optimizer=opt,metrics=['accuracy'])

    return 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 = cwd + os.path.sep + 'model'+ os.path.sep +'data'
train_path = cwd + os.path.sep + 'model'+ os.path.sep +'data'
model_path = cwd + os.path.sep + 'model'+ os.path.sep +'model_PRISMA.h5' # SPECIFICARE

# Loading the dataset
#pkl_file = open('data/data_AIRBORNE_PRISMAres.pkl', 'rb')
#pkl_file = open('data/data_AIRBORNE.pkl', 'rb')
pkl_file = open('model/data/data.pkl', 'rb')
data1 = pickle.load(pkl_file)
pkl_file.close()
data = data1[0]
X_train = data/np.max(data)
y_train = data1[1]

# Loading the dataset for the test
# = open('data/data_AIRBORNE_test_PRISMAres.pkl', 'rb')
#pkl_file = open('data/data_AIRBORNE_test.pkl', 'rb')
pkl_file = open('model/data/data_test.pkl', 'rb')
data1 = pickle.load(pkl_file)
pkl_file.close()
X_test = data1[0]
X_test = X_test/np.max(X_test)
y_test = data1[1]

print(y_train.shape,y_test.shape)

# model = load_model(model_path + 'model_keras2.h5') -> these are the best weights for Silburt

# reshaping images before feeding them to the network


train_dataset = tf.data.Dataset.from_tensor_slices((X_train,y_train))
print(train_dataset)


test_dataset = tf.data.Dataset.from_tensor_slices((X_test,y_test))
print(test_dataset)


BATCH_SIZE = 5
SHUFFLE_BUFFER_SIZE = 1000

train_batches = train_dataset.shuffle(SHUFFLE_BUFFER_SIZE).batch(BATCH_SIZE)
test_batches = test_dataset.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.15, FL=3, n_filters=60, 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]

# Initilization with previous weights
#model.load_weights('D:/Documenti_Lavoro/Craters_detection/da_francesco/moon_model_francesco.h5')

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

model = model_1D(dim=L, learn_rate=1e-4, lmbda=1e-5, drop=0.15, FL=3, n_filters=60, init='he_normal')
model.load_weights(model_path) # SPECIFICARE
eval = model.evaluate(test_batches,
                     verbose=1)
print(eval)



pred = model.predict(test_batches)
save(ref_path + 'PRISMA_CNR_predictions.npy',pred) # SPECIFICARE
print(y_test)
pred_classes = np.squeeze(np.round(pred).astype(int))
print(type(pred_classes))
print(pred_classes)
print(classification_report(y_test,pred_classes))

# pred = model.predict(test_imgs[])
# plt.imshow(pred[0],cmap='gray') -> remember to index pred


1 Physical GPUs, 1 Logical GPUs
(1905,) (272,)
<TensorSliceDataset shapes: ((234, 1), ()), types: (tf.float64, tf.int32)>
<TensorSliceDataset shapes: ((234, 1), ()), types: (tf.float64, tf.int32)>
234
Fitting the model...
Epoch 1/30
381/381 [==============================] - 5s 6ms/step - loss: 0.1610 - accuracy: 0.8068 - val_loss: 0.0796 - val_accuracy: 0.9412
Epoch 2/30
381/381 [==============================] - 1s 4ms/step - loss: 0.0353 - accuracy: 0.9877 - val_loss: 0.0153 - val_accuracy: 0.9926
Epoch 3/30
381/381 [==============================] - 1s 4ms/step - loss: 0.0130 - accuracy: 0.9969 - val_loss: 0.0079 - val_accuracy: 1.0000
Epoch 4/30
381/381 [==============================] - 1s 3ms/step - loss: 0.0065 - accuracy: 1.0000 - val_loss: 0.0105 - val_accuracy: 1.0000
Epoch 5/30
381/381 [==============================] - 1s 4ms/step - loss: 0.0048 - accuracy: 1.0000 - val_loss: 0.0346 - val_accuracy: 1.0000
Epoch 6/30
381/381 [==============================] - 1s 4ms/step - loss: 0.0050 - accuracy: 1.0000 - val_loss: 0.0124 - val_accuracy: 1.0000
Epoch 7/30
381/381 [==============================] - 1s 3ms/step - loss: 0.0034 - accuracy: 0.9989 - val_loss: 0.0617 - val_accuracy: 0.9706
Epoch 8/30
381/381 [==============================] - 1s 3ms/step - loss: 0.0026 - accuracy: 0.9996 - val_loss: 0.0088 - val_accuracy: 1.0000
Epoch 9/30
381/381 [==============================] - 1s 3ms/step - loss: 0.0025 - accuracy: 0.9999 - val_loss: 0.0391 - val_accuracy: 1.0000
Epoch 10/30
381/381 [==============================] - 1s 3ms/step - loss: 0.0026 - accuracy: 0.9995 - val_loss: 0.0997 - val_accuracy: 0.8456
Epoch 11/30
381/381 [==============================] - 1s 3ms/step - loss: 0.0033 - accuracy: 0.9976 - val_loss: 0.1325 - val_accuracy: 0.8051
Epoch 12/30
381/381 [==============================] - 1s 3ms/step - loss: 0.0018 - accuracy: 0.9998 - val_loss: 0.0875 - val_accuracy: 0.8640
Epoch 13/30
381/381 [==============================] - 1s 3ms/step - loss: 0.0048 - accuracy: 0.9949 - val_loss: 0.0251 - val_accuracy: 1.0000
Epoch 00013: early stopping
55/55 [==============================] - 0s 2ms/step - loss: 0.0066 - accuracy: 1.0000
[0.00787720363587141, 1.0]
[0 1 1 1 1 1 1 1 1 0 1 0 1 1 1 1 1 1 1 1 0 0 1 0 1 1 0 1 1 0 1 0 1 0 1 1 1
 1 0 1 1 1 1 1 1 1 0 1 1 0 0 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 0 1
 1 1 0 1 1 1 1 1 1 1 0 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 0 0 1
 0 1 1 1 1 0 1 1 1 1 1 0 1 0 1 1 1 1 0 0 1 1 1 0 1 0 0 0 1 0 1 1 1 1 0 0 1
 0 1 1 1 0 1 1 1 1 0 0 0 1 1 1 1 0 1 1 1 1 1 0 1 0 1 0 1 1 1 1 0 1 1 0 1 1
 0 0 1 1 1 1 0 1 0 1 1 1 1 1 0 0 0 0 1 1 1 1 1 1 1 0 1 1 1 1 0 1 1 0 0 1 1
 1 1 1 1 0 1 1 0 0 1 1 1 0 0 0 1 0 0 1 1 1 1 0 1 0 1 0 0 1 0 1 1 0 1 1 0 0
 1 0 0 1 0 1 0 1 1 1 1 0 0]
<class 'numpy.ndarray'>
[0 1 1 1 1 1 1 1 1 0 1 0 1 1 1 1 1 1 1 1 0 0 1 0 1 1 0 1 1 0 1 0 1 0 1 1 1
 1 0 1 1 1 1 1 1 1 0 1 1 0 0 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 0 1
 1 1 0 1 1 1 1 1 1 1 0 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 0 0 1
 0 1 1 1 1 0 1 1 1 1 1 0 1 0 1 1 1 1 0 0 1 1 1 0 1 0 0 0 1 0 1 1 1 1 0 0 1
 0 1 1 1 0 1 1 1 1 0 0 0 1 1 1 1 0 1 1 1 1 1 0 1 0 1 0 1 1 1 1 0 1 1 0 1 1
 0 0 1 1 1 1 0 1 0 1 1 1 1 1 0 0 0 0 1 1 1 1 1 1 1 0 1 1 1 1 0 1 1 0 0 1 1
 1 1 1 1 0 1 1 0 0 1 1 1 0 0 0 1 0 0 1 1 1 1 0 1 0 1 0 0 1 0 1 1 0 1 1 0 0
 1 0 0 1 0 1 0 1 1 1 1 0 0]
              precision    recall  f1-score   support

           0       1.00      1.00      1.00        81
           1       1.00      1.00      1.00       191

    accuracy                           1.00       272
   macro avg       1.00      1.00      1.00       272
weighted avg       1.00      1.00      1.00       272