Performing spatial analysis on PRISMA L2D data with prismatools¶

This notebook demonstrates how to perform spatial analysis on PRISMA hyperspectral L2D data with prismatools. For more information about PRISMA, please visit the links below:

• https://www.asi.it/en/earth-science/prisma/
• https://en.wikipedia.org/wiki/PRISMA_(spacecraft)

Note:

• Register and obtain the PRISMA datasets from the ASI PRISMA portal
• Once you have downloaded a dataset, replace the filepath which we will use in this notebook.

In [1]:

from prismatools import read_prismaL2D
from prismatools.spatial import band_math, pca_image

from prismatools import read_prismaL2D from prismatools.spatial import band_math, pca_image

In [2]:

file = "/path/to/prisma.he5"  # example: r"..\..\data\PRS_L2D_STD_20240429095823_20240429095827_0001.he5"

file = "/path/to/prisma.he5" # example: r"..\..\data\PRS_L2D_STD_20240429095823_20240429095827_0001.he5"

01. Index computation (band math)¶

Read the PRISMA dataset

In [3]:

# ds = read_prismaL2D(file_path=file)
# print(ds)

# ds = read_prismaL2D(file_path=file) # print(ds)

Compute NDVI

In [4]:

# wavelength_dict = {'red': 660, 'nir': 840}

# expressions =  {'NDVI' :'(nir - red) / (nir + red)'}

# constants = {}

# result = band_math(ds, expressions, wavelength_dict, constants, inplace=True)

# print(result)
# result.NDVI.plot(cmap='YlGn')

# wavelength_dict = {'red': 660, 'nir': 840} # expressions = {'NDVI' :'(nir - red) / (nir + red)'} # constants = {} # result = band_math(ds, expressions, wavelength_dict, constants, inplace=True) # print(result) # result.NDVI.plot(cmap='YlGn')

Compute multiple VIs

In [5]:

# wavelength_dict = {'red': 660, 'green': 560, 'blue': 480, 'nir': 840}

# expressions = {
#      'NDVI' : '(nir - red) / (nir + red)',
#      'GNDVI': '(nir - green) / (nir + green)',
#      'EVI'  : '2.5 * (nir - red) / (nir + 6 * red - 7.5 * blue + 1)', # constant might be provided directly in the expression
#      'SAVI' : '((nir - red) * (1 + L)) / (nir + red + L)'
#  }
# # or in dictionary
# constants = {
#     'L': 0.5,
#     }

# result = band_math(ds, expressions, wavelength_dict, constants, inplace=True)
# print(result)

# wavelength_dict = {'red': 660, 'green': 560, 'blue': 480, 'nir': 840} # expressions = { # 'NDVI' : '(nir - red) / (nir + red)', # 'GNDVI': '(nir - green) / (nir + green)', # 'EVI' : '2.5 * (nir - red) / (nir + 6 * red - 7.5 * blue + 1)', # constant might be provided directly in the expression # 'SAVI' : '((nir - red) * (1 + L)) / (nir + red + L)' # } # # or in dictionary # constants = { # 'L': 0.5, # } # result = band_math(ds, expressions, wavelength_dict, constants, inplace=True) # print(result)

In [6]:

# # plotting
# result["SAVI"].plot(cmap='YlGn', vmin=0, vmax=1)

# # exporting
# indices = [v for v in result.data_vars if v != 'reflectance']
# for index in indices:
#     output_path = fr"C:\Users\loren\Desktop\prismatools\out_test\{index}.tif"
#     result[index].rio.to_raster(output_path)

# # plotting # result["SAVI"].plot(cmap='YlGn', vmin=0, vmax=1) # # exporting # indices = [v for v in result.data_vars if v != 'reflectance'] # for index in indices: # output_path = fr"C:\Users\loren\Desktop\prismatools\out_test\{index}.tif" # result[index].rio.to_raster(output_path)

02. Performing Principal Component Analysis (PCA)¶

• From GeoTiff file (from path)

In [7]:

input_file = "path/to/TIFF_image.tif"  # replace with actual path
output_file = "path/to/out_TIFF_image.tif"  # replace with actual path

input_file = "path/to/TIFF_image.tif" # replace with actual path output_file = "path/to/out_TIFF_image.tif" # replace with actual path

In [8]:

# n_components = 10

# ds_pca, expl_var = pca_image(input_file, output_file, n_components=n_components)
# print(ds_pca)
# print(expl_var)

# n_components = 10 # ds_pca, expl_var = pca_image(input_file, output_file, n_components=n_components) # print(ds_pca) # print(expl_var)

In [9]:

# ds_pca["pca"].sel(pc=2).plot(x="x", y="y")

# ds_pca["pca"].sel(pc=2).plot(x="x", y="y")

• From xarray.Dataset

In [10]:

# n_components = 10

# ds_pca, expl_var = pca_image(ds, output_file, n_components=n_components)
# print(ds_pca)
# print(expl_var)

# n_components = 10 # ds_pca, expl_var = pca_image(ds, output_file, n_components=n_components) # print(ds_pca) # print(expl_var)