Pansharpening Landsat con Qgis

 

Si caricano i raster delle bande 2,3,4 (BGR) e la banda 8 pancromatica

Si apre il menu Raster/Miscellanea/Crea Raster Virtuale selezionando come input le bande ottiche e spuntando il flag Place each input file in a separate band

A questo punto dal toolbox si seleziona pansharpening indicando come Dataset spettrale il raster virtuale creato in precedenza e come dataset pancromatico la banda 8

Snap e Landsat

 Per utilizzare le immagini Landsat all'interno del software Snap di Esa e' sufficiente spacchettare il file tar ed aprire il file txt con contiene nel nome la sigla MTL

RGB change detection con RasterIO

 Un semplice metodo per avere il change detection da una immagine RGB e' la formula seguente

 

 che non e' altro che la distanza euclidea tra i due punti rappresentativi delle immagini nello spazio delle bande (quindi puo' essere scalato anche su n-bande)

 

 utilizzando la libreria RasterIO si puo' scrivere il codice seguente

 

import numpy as np

import rasterio

from rasterio.plot import show_hist,show

immagine1 = "cd1.tif"

immagine2 = "cd2.tif"

im1 = rasterio.open(immagine1)

im2 = rasterio.open(immagine2)

im1_b2 = im1.read(1).astype(float)

im1_b3 = im1.read(2).astype(float)

im1_b4 = im1.read(3).astype(float)

im2_b2 = im2.read(1).astype(float)

im2_b3 = im2.read(2).astype(float)

im2_b4 = im2.read(3).astype(float)

np.seterr(divide = "ignore", invalid = "ignore")

change_det = np.zeros(im1_b2.shape, dtype=rasterio.float32)

change_det = np.sqrt(np.power(im1_b2-im2_b2,2)+np.power(im1_b3-im2_b3,2)+np.power(im1_b4-im2_b4,2))

#print("Change Det Max " + str(np.max(change_det)))

#print("Change Det Min " + str(np.min(change_det)))

#show(change_det)

#show_hist(change_det)

kwargs = im1.meta

kwargs.update(driver='GTiff',dtype=rasterio.float32,count=1)

cd_image = rasterio.open('CD.tiff','w',**kwargs )

cd_image.write(change_det,1)

cd_image.close()

 

L'importanza di una antenna UV-5118

 Mi sono comprato su Aliexpress questa radio UV-5118 che era venduta con Air Band...una volta arrivata non solo non riceveva in air band ma praticamente a nessuna frequenza 

Dopo aver ottenuto il rimborso ho sostituito l'antenna con una presa da un Baofeng...la air band non si riceve ma in UHF e VHF ha iniziato a funzionare in modo decisamente buono

Satellite Derived Bathymetry Mapping

Una applicazione dello script a questo link applicato direttamente in EO Browser

Versilia

Booca L'Arno

Lago di massaciuccoli

Golfo di Spezia Foce del Magra

in linea di principio lo script puo' essere semplificato dato che la batimetria puo' essere stimata in modo relativa con il rapporto log10(B3/B4) di Sentinel 2 (attenzione...Rrs non e' la Surface Reflectance che si trova nel prodotto L2 di Sentinel ma e' la Remote Sensing Reflectance)

in ogni caso si vede che l'influenza del carico di sedimenti influenza la stima della batimetria (vedi immagine successiva. Batimetria ripresa da progetto Camp Italy)

Correlazione tra Sentinel 5P e stazioni a terra

In questo post, usando gli script di Earth Engine visti in precedenza, ho provato a mettere in correlazione il dato telerilevato di Sentinel 5P con alcune stazioni di monitoraggio aria Arpat (dati pubblici sul sito www.arpat.toscana.it) 

Parametro NO2

Vi e' un ottimo accordo tra i dati medi mensili delle concentrazioni di NO2 da satellite ed a terra

la curva di correlazione non e' incoraggiante

Parametro CH4

I dati da satellite mostrano il gradiente di aumento del metano che e' coerente con il valore di crescita mondiale 

Parametro S02

Per il parametro S02 non si possono fare correlazioni perche' quando le concentrazioni aumentano (nel periodo autunnale) il satellite non riporta dati

Stazione Arpat via Bassi Firenze

Mappa SO2 Piana Fiorentina Giugno 2019

Sentinel 5P

Parametro O3

Per l'ozono ho molti dubbi...comunque esporti i dati si vede che c'e' un offset di circa 3 mesi tra i dati telerilevati e quelli a terra

Parametro CO

Questo parametro e' critico per la correlazione con la stazione a terra perche' a Firenze la stazione di misura non e' in posizione di fondo ma a bordo di una viabilita' principale

Parametro PM10

Il parametro PM10 puo' essere derivato da UV Aerosol nel dataset di Sentinel 5P

Qui il  confronto delle medie mensili tra i dati satellitari ed una stazione a terra

Aerosol marzo 2022

qui un riferimento bibliografico della correlazione tra dati satellitari ed a terra

Earth Engine Sentinel 5P time series NO2

 Due script per scaricare serie tempo dei dati Sentinel 5P da Google Earth Engine dato un punto geografico

Tutti i dati

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

var point = ee.Geometry.Point([11.097849, 43.793234]);
var sentinel = ee.ImageCollection('COPERNICUS/S5P/NRTI/L3_NO2');
var sentinelLST = sentinel.filterBounds(point)
.filterDate('2019-01-01', '2022-12-31')
.select('NO2_column_number_density');

sentinelLST = sentinelLST.map(function(img){
var date = img.get('system:time_start');
return img.multiply(100000).set('system_time_start', date);
});

var createTS = function(img){
var date = img.get('system_time_start');
var value = img.reduceRegion(ee.Reducer.mean(), point).get('NO2_column_number_density');
var ft = ee.Feature(null, {'system:time_start': date,
'date': ee.Date(date).format('Y/M/d'),
'value': value});
return ft;
};

var TS = sentinelLST.map(createTS);

var graph = ui.Chart.feature.byFeature(TS, 'system:time_start', 'value');

print(graph.setChartType("Sentinel 5P")
.setOptions({vAxis: {title: 'NO2'},
hAxis: {title: 'Date'}}));

Export.table.toDrive({collection: TS, selectors: 'date, value'});

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

Medie Giornaliere

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

var collection = ee.ImageCollection('COPERNICUS/S5P/OFFL/L3_NO2')
  .select('tropospheric_NO2_column_number_density')
var daily_data = ee.ImageCollection(ee.List.sequence(2019,2019).map(function(year){
  var date1 = ee.Date.fromYMD(year,1,1)
  var date2 = date1.advance(1,'year')
  //Calculate the number of days per year
  var doy = date2.difference(date1,'day')
  var doyList = ee.List.sequence(1,doy)
  //Daily image mean synthesis using doy
  var day_imgs = doyList.map(function(doy){
    doy = ee.Number(doy)
    var temp_date = date1.advance(doy.subtract(1),"day")
    var temp_img = collection.filterDate(temp_date,temp_date.advance(1,'day'))
    return temp_img.mean().set("system:time_start",temp_date.millis())
  })
  return day_imgs
}).flatten())
Map.addLayer(daily_data)
var chart = ui.Chart.image.series({
imageCollection:daily_data,
region:roi,
reducer:ee.Reducer.mean(),
scale:1113.2,
// xProperty:,
})
print(chart)

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

Medie Mensili

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

var point = ee.Geometry.Point([11.097849, 43.793234]);
var dataset = ee.ImageCollection('COPERNICUS/S5P/NRTI/L3_NO2')
          .filterBounds(point)
          .filterDate('2019-01-01', '2022-12-31')
          .select('NO2_column_number_density');
var months = ee.List.sequence(1, 12);
var start_year = 2019;
var start_date = '2019-01-01';
var end_year = 2022;
var end_date = '2022-12-31';
var years = ee.List.sequence( start_year, end_year);

var byMonthYear =  ee.FeatureCollection(
  years.map(function (y) {
    return months.map(function(m){
          var w = dataset.filter(ee.Filter.calendarRange(y, y, 'year'))
                    .filter(ee.Filter.calendarRange(m, m, 'month'))
                    .mean();
      var pointMean = w.reduceRegion({reducer:ee.Reducer.first(), geometry:point,scale:1000});  
      return ee.Feature(null).set("value",pointMean.get("NO2_column_number_density")).set("year",y).set("month",m);
    })
  }).flatten()
);

print("feature collection",byMonthYear);

Export.table.toDrive({collection:byMonthYear,description:"O3_signa"})

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