Primi passi con R PCA

Si parte da un file CSV in cui le colonne sono le variabili (in questo caso dati analitici) 

La prima riga e' di intestazione, separatore punto e virgola, punto decimake

====================================================

Punto;As;Co;Hg
Pz3A;26.4;13.4;0.206
Pz3B;41.9;17.7;0.281
Pz4;8.2;6.05;1.13
Pz5A;21.5;18.6;0.281
Pz5B;20.1;18.2;0.235
Pz6;5.8;20.4;0.138
Pz7A;6.8;23.2;0.082
Pz7B;7.4;24;0.090
Pz8A;16.8;21.8;0.325
Pz8B;12.4;22.2;0.279
Pz9A;46.9;30.2;0.71
Pz9B;20.0;21.6;0.319
Pz10A;29.3;9.4;1.03
Pz10B;18.5;56.1;0.45
Pz11A;63;12.7;0.210
Pz12;32.1;3.53;0.0298
Pz13A;18.2;31.1;0.75
Pz13B;14.8;29.1;0.200
Pz14A;9.5;28.6;0.043
Pz14B;9.2;32.6;0.045

====================================================

Per elaborare i dati con PCA vengono caricate le librerie FactoMineR e factoextra

Per leggere il file CSV si deve specificare il separatore decimale

Si deve anche indicare che la prima colonna e' un dato qualitativo (da escludere dal calcolo), le dimensioni della PCA pari a 3,

 

====================================================

install.packages("FactoMineR")
install.packages("factoextra")

library("FactoMineR")
library("factoextra")
data_res <- read.csv2("c:/Users/l.innocenti/Desktop/arsenico/arsenico.csv",dec=".")
data_res.pca <- PCA(data,scale.unit = TRUE,ncp=3,graph = TRUE,quali.sup=1)
data_res.pca
eig.val <- get_eigenvalue(data_res.pca)
fviz_eig(data_res.pca,addlabels = TRUE)
fviz_pca_var(data_res.pca)
var <-get_pca_var(data_res.pca)
library(corrplot)
corrplot(var$cos2,is.corr = FALSE)
fviz_pca_ind(data_res.pca)

Eliminare dati a bassa coerenza da inteferogrammi Sentinel 1

Nel caso i dati di coerenza di un inteferogramma sentinel 1 siano bassi (tali da rendere inutili ulteriori elaborazioni) si puo' tentare di filtrare i dati di fase con bassa coerenza

 

Per fare cio' si seleziona la banda della fase e con tasto destro si selezione Properties. In seguito si edita Pixel Values Expression con una formula tipo quella sottostante

if coh_IW1_VV_12Dec2022_24Dec2022 > 0.4 then atan2(atan2(q_ifg_IW1_VV_12Dec2022_24Dec2022,i_ifg_IW1_VV_12Dec2022_24Dec2022))  else 0

in pratica tutti i pixel con coerenza inferiore a 0.4 vengono messi a zero

Si salva il prodotto e si puo' procedere con i passi successivi

DEM e coerenza su Sentinel 1

 Dopo un po' di insuccessi nella creazione di DEM da dati Sentinel 1 ho ricercato in bibliografia le possibili cause

Per prima cosa dopo aver fatto l'interferogramma si deve selezionare l'immagine di coerenza e creare l'istogramma dal menu Analysis/Histogram (premendo refresh),

In generale per avere un DEM significativo il valore di coerenza di almeno il 70% dei punti deve essere almeno 0.4

Scarsi valori di coerenza sono imputabili a cause atmosferiche

Viene indicato da vari autori un offset tra 20 e 40 m nei dati di elevazione tra il dato reale ed il dato radar

Frecce tricolori 12 febbrario 2023

 

Frecce tricolori Firenze 11 febbraio 2023

 

zona Firenze Sud

Sentinel 5 NetCDF CH4

Ho provato a scaricare i dati del Sentinel 5, un satellite specializzato nel rilevamento di gas atmosferici (CH4, CO, formaldeide, O3,...)

I dati si possono scaricare da SciHub (https://s5phub.copernicus.eu/dhus/#/home) e sono distribuiti in formato NetCDF. Le dimensione del pixel sono di circa 4x7 Km

il dato e' da intendersi come in attesa di validazione (Sentinel 5 e' in fase pre-operativa). Come indicato anche nelle note di Google Earth Engine non solo basta scartare i pixel con il valore di qa <0.5 (e sono tanti in ogni immagine). Il valore di q5 si estrae direttamente dal file netcdf

Timelapse Luglio 2022 CH4 

I dati possono essere visualizzati con Planoply ma ho preferito provare ad utilizzare le librerie di Python 

from netCDF4 import Dataset

import numpy as np

from mpl_toolkits.basemap import Basemap

import sys

file = sys.argv[1]

data = Dataset(file, mode='r') # read the data

print(type(data)) # print the type of the data

print(data.groups['PRODUCT']) # print the variables in the data

print (data.groups['PRODUCT'].variables.keys())

ch4 = data.groups['PRODUCT'].variables['methane_mixing_ratio'][0,:,:]

lons = data.groups['PRODUCT'].variables['longitude'][:][0,:,:]

lats = data.groups['PRODUCT'].variables['latitude'][:][0,:,:]

print (lons.shape)

print (lats.shape)

print (ch4.shape)

ch4_units = data.groups['PRODUCT'].variables['methane_mixing_ratio'].units

import matplotlib.pyplot as plt

from matplotlib.colors import LogNorm

from mpl_toolkits.basemap import Basemap

#lon_0 = lons.mean()

#lat_0 = lats.mean()

lon_0 = 11

lat_0 = 43

m = Basemap(width=1000000,height=1500000,

resolution='l',projection='stere',\

lat_ts=2,lat_0=lat_0,lon_0=lon_0)

xi, yi = m(lons, lats)

# Plot Data

cs = m.pcolor(xi,yi,np.squeeze(ch4),norm=LogNorm(), cmap='jet')

# Add Grid Lines

m.drawparallels(np.arange(-80., 81., 10.), labels=[1,0,0,0], fontsize=10)

m.drawmeridians(np.arange(-180., 181., 10.), labels=[0,0,0,1], fontsize=10)

# Add Coastlines, States, and Country Boundaries

m.drawcoastlines()

m.drawstates()

m.drawcountries()

# Add Colorbar

cbar = m.colorbar(cs, location='bottom', pad="10%")

cbar.set_label(ch4_units)

# Add Title

plt.title(file)

plt.get_current_fig_manager().full_screen_toggle() # toggle fullscreen mode

#plt.show()

plt.savefig(file+".png",dpi=300)

una elaborazione puo' essere fatta anche Google Earth Engine

=========================

var collection = ee.ImageCollection('COPERNICUS/S5P/OFFL/L3_CH4')

  .select('CH4_column_volume_mixing_ratio_dry_air')

  .filterDate('2019-01-01', '2019-12-31');

  

var band_viz = {

  min: 1000,

  max: 2000,

  palette: ['black', 'blue', 'purple', 'cyan', 'green', 'yellow', 'red']

};

var median = collection.reduce(ee.Reducer.median());

print(median)

Map.addLayer(median, band_viz, 'S5P CH4');

Map.setCenter(11.295553516712024, 43.976202446636606,10);

PCA e GLCM in Sentinel 1

Utilizzando le immagini GRDH di Sentinel 1 (il dato GRDH si differenzia dal SLC perche' nel primo si ha solo l'ampiezza e l'intesita' del segnale nelle due polarizzazioni VV e VH mentre il SLC si ha sia la componente reale che quella complessa) si possono utilizzare i livelli di risposta (scale di grigi nelle immagini) per applicare degli algoritmi di PCA (componenti principali) e GLCM (Gray Level Co occurence matrix) per analisi delle texture

Dissimilarity map da GLCM

Composit di componenti principali

Due esempi di elaborazione in PCA e GLCM pre e post evento di inondazione

Partendo dal dato GRDH i passi da seguire sono (https://www.youtube.com/watch?v=xoQ4NikdOq0)

• Subset
• Apply Orbit File
• Radiometric -> Calibrate
• Speckle Filtering
• Range Doppler Terrain Correction
• Raster/Data Conversion/Convert to dB
• Raster/Image Analysis/ Principal components
• Raster/Image Analysis/Texture Analysis/Gray Level Co occurence Matrix 

Nel caso di PCA nelle bande troveremo le mappe per ogni componente principale con in aggiunta la mappa di errore e  response

Nel caso di GCLM troveremo per sigma0 VV e VH la mappa di contrasto, dissimilarita', omogeneita', ASM, energia, entropia, media, varianza,correlazione