Visualizzazione post con etichetta Google Earth Engine. Mostra tutti i post
Visualizzazione post con etichetta Google Earth Engine. Mostra tutti i post

venerdì 21 aprile 2023

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"})

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


mercoledì 8 marzo 2023

Formaldeide e metano da Sentinel 5 con Google Earth Engine

 Usare Google Earth Engine per Sentinel 5 e' comodo perche' i pixel sono stati ricalcolati in modo omogeneo (nei file netcdf originali i pixel hanno dimensioni differenti a secondo dell'orbita)

I valori di CH4 sono espressi in ppm ed il range di valore e' estremamente stretto (da 1860 a 1890 ppm nell'immagine ) con la dimensione del pixel di 1113 m. I dati sono disponibili sono per le aree di fondovalle




I dati di formaldeide sono espressi in mol/m2 ed il range di valori varia tra 0.00007 e 0.00014 con pixel di 1113 m. 




var palettes = require('users/gena/packages:palettes');
  
var palette_r = palettes.colorbrewer.Accent[5];
var collection = ee.ImageCollection('COPERNICUS/S5P/NRTI/L3_HCHO')
  .select('tropospheric_HCHO_column_number_density')
  .filterDate('2019-01-01', '2019-12-31');
var median = collection.reduce(ee.Reducer.mean());
var band_viz = {
  min: 0.00007,
  max: 0.00014,
  palette: ['black', 'blue', 'purple', 'cyan', 'green', 'yellow', 'red']
};
Map.addLayer(collection.mean(), band_viz, 'S5P HCHO');
Map.setCenter(11.295553516712024, 43.976202446636606,9);
Map.style().set('cursor', 'crosshair');
var header = ui.Label('HCHO', {fontWeight: 'bold', fontSize: '18px'});
var toolPanel = ui.Panel([header], 'flow', {width: '400px'});
Map.onClick(function(coords) {
  // Create or update the location label (the second widget in the panel)
  var location = 'lon: ' + coords.lon.toFixed(4) + ' ' +
                 'lat: ' + coords.lat.toFixed(4);
  var click_point = ee.Geometry.Point(coords.lon, coords.lat);
  var CH4Value = median.reduceRegion(ee.Reducer.first(), click_point, 90).evaluate(function(val){
    console.log(val);
    var CH4Text = 'CH4: ' + val;
    toolPanel.widgets().set(2, ui.Label(CH4Text));
  });
  
  toolPanel.widgets().set(1, ui.Label(location));
  // Add a red dot to the map where the user clicked.
  Map.layers().set(1, ui.Map.Layer(click_point, {color: 'FF0000'}));
});
// Add the panel to the ui.root.
ui.root.add(toolPanel);

giovedì 17 novembre 2022

Geotiff da assets in EarthEngine

In Earth Engine e' possibile visualizzare dei Geotiff caricandoli in assets 


I geotiff devono essere puri e non tiff con allegato il world file in formato tfw.

Per la conversione da tfw a geotiff si puo' usare Gdal Translate includendo anche il sistema di riferimento

(quello che segue e' un esempio di conversione di un foglio CARG Carta geologica regionale Toscana)

gdal_translate -a_srs EPSG:3003 -of GTiff 264020.tif 264020_.tif

Le immagini che sono in assets non sono immediatamente leggibili dalla versione app del progetto.

Devono essere esplicitamente impostati i permessi di lettura nelle proprieta' dell'asset


Per visualizzare file prendendolo dagli assets si usa la seguente sintassi

var geologia = ee.Image('projects/ee-lucainnoc/assets/264020');

Map.addLayer(geologia, {},'Geologia',false,0.5);


mercoledì 16 novembre 2022

Soil Moisture Sentinel 1 Google Earth Engine

 Visto che le immagini radar presenti su Google EarthEngine sono gia' radiometricamente e geometricamente corrette ho voluto provare a replicare su questa piattaforma quanto visto nel precedente post



la demo funzionante e' all'indirizzo sotto riportato

https://lucainnoc.users.earthengine.app/view/soil-moisture-sentinel-1

Questa e' la spvrapposizione tra la geologia e l'algoritmo dell'area del precedente post (Casole, Vicchio di Mugello). 

La mappa e' la differenza di umidita' tra il periodo primaverile ed estivo dell'anno 2021 con il colore blu che indica massima variazione negativa e il colore rosso massima variazione positiva


function maskS2clouds(image) {
  var qa = image.select('QA60');
  // Bits 10 and 11 are clouds and cirrus, respectively.
  var cloudBitMask = 1 << 10;
  var cirrusBitMask = 1 << 11;
  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0)
      .and(qa.bitwiseAnd(cirrusBitMask).eq(0));
  return image.updateMask(mask).divide(10000);
}

var imgVV = ee.ImageCollection('COPERNICUS/S1_GRD')
        .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))
        .filter(ee.Filter.eq('instrumentMode', 'IW'))
        .select('VV')
        .map(function(image) {
          var edge = image.lt(-30.0);
          var maskedImage = image.mask().and(edge.not());
          return image.updateMask(maskedImage);
        });
//seleziona solo ASCENDENTE
var asc = imgVV.filter(ee.Filter.eq('orbitProperties_pass', 'ASCENDING'));

//seleziona il periodo primaverile ed estivo
var spring = ee.Filter.date('2021-03-01', '2021-04-20');
var summer = ee.Filter.date('2021-06-11', '2021-08-31');
//seleziona il periodo dal gennaio 2019 al dicembre 2021
//per ottenere i valori massimi e minimi
var full = ee.Filter.date('2019-01-01', '2021-12-31');
var sar_dataset_min = ee.Image.cat(asc.filter(full).min());
var sar_dataset_max = ee.Image.cat(asc.filter(full).max());
// calcola la variazione massima nel periodo di interesse 2019-2021
var sensitivity = sar_dataset_max.subtract(sar_dataset_min);

var primavera = ee.Image.cat(asc.filter(spring).mean());
var estate = ee.Image.cat(asc.filter(summer).mean());
var Mv_primavera = (primavera.subtract(sar_dataset_min)).divide(sensitivity);
var Mv_estate = (estate.subtract(sar_dataset_min)).divide(sensitivity);
var delta = Mv_estate.subtract(Mv_primavera);
//var Mv_estate = (estate - sar_dataset_min)/(sensitivity);
//print(Mv_primavera);
//estre i truecolor ottici per riferimento
var rgb_spring = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')
                  .filterDate('2018-03-01', '2018-04-20')
                  // Pre-filter to get less cloudy granules.
                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE',20))
                  .map(maskS2clouds);
var visualization = {
  min: 0.0,
  max: 0.3,
  bands: ['B4', 'B3', 'B2'],
};
var rgb_summer = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')
                  .filterDate('2018-06-11', '2018-08-31')
                  // Pre-filter to get less cloudy granules.
                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE',20))
                  .map(maskS2clouds);
var visualization = {
  min: 0.0,
  max: 0.3,
  bands: ['B4', 'B3', 'B2'],
};

Map.setCenter(11.489351719442514,43.95792446517467, 14);
//Map.addLayer(primavera, {min: -25, max: 5}, 'Spring', true);
//Map.addLayer(estate, {min: -25, max: 5}, 'Summer', true);
Map.addLayer(rgb_spring.mean(), visualization, 'RGB Spring');
Map.addLayer(rgb_summer.mean(), visualization, 'RGB Summer',false);
var SensitivityBandVis = {
  min: 0,
  max: 50,
  palette: ['blue', 'yellow', 'green']
};
//Map.addLayer(rgb_spring.mean(), visualization, 'RGB Spring');
//Map.addLayer(sensitivity, SensitivityBandVis, 'Sensitivity');

var MvBandVis = {
  min: 0,
  max: 1,
  palette: ['blue', 'yellow', 'green']
};

var DeltaBandVis = {
  min: -1,
  max: 1,
  palette: ['blue', 'yellow', 'green']
};
Map.addLayer(Mv_primavera, MvBandVis, 'Mv Spring',false);
Map.addLayer(Mv_estate, MvBandVis, 'Mv Summer',false);
Map.addLayer(delta, DeltaBandVis, 'Mv Delta',true,0.5);

martedì 23 agosto 2022

Copertura nevosa con Google Earth Engine ed ERA5

Questo dato e' stato ricercato perche' in un settore dell'Appennino ToscoEmiliano qualche anno fa si sono riattivate numerose frane. Il sospetto e' che fossero cadute nevicate tardive che si sono poi sciolte rapidamente con l'arrivo del caldo invernale innescando i movimenti di versante

Sull'area non sono presenti nivometri per cui non vi sono dati di verita' a terra. Le immagini Landsat ovviamente nel periodo invernale sono spesso nuvolose e di poco utilizzo

Un primo metodo e' quello di utilizzare i dati di ERA5

Questa analisi si puo' fare interrogando direttamente le API di Copernicus senza passare da Google Earth Engine

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

var geometry = /* color: #d63000 */ee.Geometry.MultiPoint(
        [[10.2014189280409, 44.47706314233407],
         [10.389914585837579, 44.524348579528535],
         [10.228782890268272, 44.53074818785566],
         [10.075498095653819, 44.458900035551736],
         [10.100715972012257, 44.51255997214677]]);
/***** End of imports. If edited, may not auto-convert in the playground. *****/

var point = ee.Geometry.Point([10.2014189280409, 44.47706314233407]);
var snow_depth = ee.ImageCollection('ECMWF/ERA5_LAND/MONTHLY')
            .select('snow_depth')
            .filterBounds(point)
            .filter(ee.Filter.calendarRange(2000,2020,'year'))
            .map(function(image){return image.clip(point)}) ;
      

// plot full time series
print(
  ui.Chart.image.series({
    imageCollection: snow_depth,
    region: point,
    scale: 1000
  }).setOptions({title: 'Snow Monthly ERA5'})
);

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

Altrimenti si puo' usare l'indice NDIS_Snow_Cover di Modis che ha un dettagli temporale molto maggiore


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

var geometry = /* color: #d63000 */ee.Geometry.MultiPoint(
        [[10.2014189280409, 44.47706314233407],
         [10.389914585837579, 44.524348579528535],
         [10.228782890268272, 44.53074818785566],
         [10.075498095653819, 44.458900035551736],
         [10.100715972012257, 44.51255997214677]]);
/***** End of imports. If edited, may not auto-convert in the playground. *****/


var startDate = ee.Date('2010-01-01'); // set start time for analysis
var endDate = ee.Date('2022-01-01'); // set end time for analysis
var nDays = ee.Number(endDate.difference(startDate,'day')).round();

var point = ee.Geometry.Point([10.04337, 45.71931]);

var chirps = ee.ImageCollection('MODIS/006/MOD10A1')
            .select('NDSI_Snow_Cover')
            .filterBounds(point)
            .filterDate(startDate, endDate)
            .map(function(image){return image.clip(point)}) ;

      



var byday = ee.ImageCollection(
  // map over each day
  ee.List.sequence(0,nDays).map(function (n) {
    var ini = startDate.advance(n,'day');
    // advance just one day
    var end = ini.advance(1,'day');
    return chirps.filterDate(ini,end)
                .select(0).mean()
                .set('system:time_start', ini);
}));


//print(byday);
Map.setCenter(10.04582, 45.71337, 11);

//Map.addLayer(ee.Image(byday.mean()),{min: 0, max: 50},'Precipitazioni');

// plot full time series
print(
  ui.Chart.image.series({
    imageCollection: byday,
    region: point,
    scale: 1000
  }).setOptions({title: 'MODIS NDSI'})
  );

Landsat 8 truecolor gif con Google Earth Engine

 Una prova di LandSat 8 truecolor con Google Earth Engine...stavo cercando i periodi con neve su questo settore dell'Appennino ma non e' mi e' stato molto di aiuto (il periodo va da agosto 2013 a agosto 2022....28 immagini con lista completa al termine del post)


var geometry = 
    /* color: #d63000 */
    /* displayProperties: [
      {
        "type": "rectangle"
      }
    ] */
    ee.Geometry.Polygon(
        [[[10.11590875374765, 44.54755512745409],
          [10.11590875374765, 44.39388022912829],
          [10.329541337365814, 44.39388022912829],
          [10.329541337365814, 44.54755512745409]]], null, false);
var visualization = {"bands":["B4","B3","B2"],"min":0,"max":0.3};
// Applies scaling factors.
function applyScaleFactors(image) {
  var opticalBands = image.select('SR_B.').multiply(0.0000275).add(-0.2);
  var thermalBands = image.select('ST_B.*').multiply(0.00341802).add(149.0);
  return image.addBands(opticalBands, null, true)
              .addBands(thermalBands, null, true);
}
// Landat 8 surface reflection data
var L8Coll = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
    .filterBounds(geometry).sort('DATE_ACQUIRED')
    .map(applyScaleFactors)
    .filterMetadata('CLOUD_COVER', 'less_than', 1)
    .map(function(image){return image.clip(geometry)});
    
print(L8Coll)
var L8Coll=L8Coll.select(['SR_B2', 'SR_B3', 'SR_B4','SR_B5', 'SR_B6', 'SR_B7','ST_B10'] , ['B2', 'B3', 'B4','B5', 'B6', 'B7','B10'])
print(L8Coll)
Map.addLayer(L8Coll.first(), visualization, 'L8 2013 ');   
Map.addLayer(L8Coll.sort('CLOUD_COVER').first(), visualization, 'L8 Best ');   
Map.centerObject(L8Coll)
print('L8Coll Least Cloud Cover',L8Coll.sort('CLOUD_COVER').first())
print('Landsat Collection Date Acquired',L8Coll.aggregate_array('DATE_ACQUIRED'))
print('Landsat Collection Cloud Cover',L8Coll.aggregate_array('CLOUD_COVER'))
print('Landsat Collection Path/Row',L8Coll.aggregate_array('WRS_PATH'),
      L8Coll.aggregate_array('WRS_ROW'))
var L8LeastCC=L8Coll.sort('CLOUD_COVER').first()
print('L8 Least Cloud Cover',L8LeastCC)
print(ee.Date(L8Coll.first().get('system:time_start')).format("yyyy-MM-dd"))
//+++++++++++++++++++++
/*
// Load a Landsat 8 collection for a single path-row.
var collection = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
  .filter(ee.Filter.eq('WRS_PATH', 179))
  .filter(ee.Filter.eq('WRS_ROW', 34));
*/
// This function adds a band representing the image timestamp.
var addTime = function(image) {
  return image.addBands(image.metadata('system:time_start'));
};
// Map the function over the collection and display the result.
var L8CollT=L8Coll.map(addTime)
print('L8coll system Time',L8CollT);

print(ee.Date(L8Coll.first().get('system:time_start')).format("yyyy-MM-dd"))
var doy = ee.Date(L8Coll.first().get('system:time_start')).getRelative('day', 'year');
print('doy',doy)
//++++++++++++++++++++++

// Create RGB visualization images for use as animation frames.
var rgbVis = L8Coll.map(function(img) {
  return img.visualize(visualization).clip(geometry)
});
print('rgbVis',rgbVis)
// Define GIF visualization parameters.
var gifParams = {
  'region': geometry,
  'dimensions': 800,
  'crs': 'EPSG:3857',
  'framesPerSecond': 1//,  'format': 'gif'
};
// Print the GIF URL to the console.
print(rgbVis.getVideoThumbURL(gifParams));
// Render the GIF animation in the console.
print(ui.Thumbnail(rgbVis, gifParams));
////+++++++++++++++

0:
2013-08-03
1:
2013-08-19
2:
2013-09-04
3:
2014-06-10
4:
2016-07-17
5:
2016-08-27
6:
2017-03-14
7:
2017-04-08
8:
2017-06-11
9:
2017-08-05
10:
2017-08-30
11:
2018-06-30
12:
2019-01-15
13:
2019-06-24
14:
2019-07-26
15:
2019-09-12
16:
2019-09-21
17:
2020-04-07
18:
2020-04-23
19:
2020-07-21
20:
2020-08-22
21:
2021-03-02
22:
2021-10-28
23:
2022-02-08
24:
2022-03-28
25:
2022-04-29
26:
2022-07-02
27:
2022-08-03

giovedì 21 luglio 2022

Google Earth Cloud Engine NOAA GFS

 Invece di procedere in modo tradizionale come indicato nel post precedente si possono interrogare i dati GFS NOAA anche da Google Earth Engine con il vantaggi che i dati sono gia' scaricati e pronti all'elaborazione (con lo svantaggio che i risultati possono essere esportati in automatico solo su cloud e non sulla macchina locale)


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

var palettes = require('users/gena/packages:palettes');
var palette = palettes.colorbrewer.RdYlGn[9];

var region = ee.Geometry(Map.getBounds(true))
// Strip time off current date/time
var today = ee.Date(new Date().toISOString().split('T')[0])
var image;
var collection = ee.ImageCollection('NOAA/GFS0P25')
  .select('temperature_2m_above_ground')
  .filterDate(today.advance(-1, 'day'), today.advance(1, 'day'))
  .sort('system:time_start', false)
console.log(collection);
// per i successivi 3 giorni le previsioni sono orarie
// le previsioni vengono ricalcolate ogni 6 ore
var forecasts = ee.ImageCollection(
  [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23, 24,
     25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,
     49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,63,64,65,66,67,68,69,70,71].map(forecast) 
  //[12, 24, 36,48,60, 72].map(forecast) 
)
Map.setCenter(11.8537, 46.4297, 12);
Map.setOptions('SATELLITE');
Map.setControlVisibility(false, false, true, true, true, false, false)
Map.addLayer(forecasts, {min: 0, max: 30,palette: palette,opacity: 0.5}, 'Temperature forecast', true)
var textStyle = {
  color: 'grey',
  fontName: 'arial',
  fontSize: 10,
  bold: true,
  italic: false
}
var chart = ui.Chart.image.series({
  imageCollection: forecasts, 
  region: geometry,
  reducer: ee.Reducer.mean(), 
  scale: forecasts.first().projection().nominalScale(),
  xProperty: 'date',
}).setChartType('LineChart').setOptions({title: 'Forecast 72h Temperature 2m above ground (C)',hAxis: {title: 'Date',textStyle: textStyle},vAxis: {title: 'Temperature'}});

//print(chart)
var panel_chart = ui.Panel({style: {width: '400px'}})
     .add(ui.Label('NOAA/GFS0P25'))
     .add(chart);
ui.root.add(panel_chart);


function forecast(hours) {
  image = collection
    .filterMetadata('forecast_hours', 'equals', hours)
    // Since colleciton is sorted descending, if there are more than 
    // one forecast, first() will give the latest.
    .first()
  var date = image.date().advance(hours, 'hours')
  return image
    .set('date', date.format())
}
var panel = ui.Panel();
panel.style().set({
  width: '400px',
  position: 'bottom-right'
});
//Map.add(panel);
Map.onClick(function(coords) {
  //panel_chart.clear();
  console.log(collection);
  //Map.add(panel);
  var location = 'lon: ' + coords.lon.toFixed(4) + ' ' +
                 'lat: ' + coords.lat.toFixed(4);
  var click_point = ee.Geometry.Point(coords.lon, coords.lat);
  var valore = image.reduceRegion(ee.Reducer.first(), click_point, 90).evaluate(function(val){
    //console.log(valore);
    //var demText = 'Habitat suitability: ' + val.b1;
    //toolPanel.widgets().set(2, ui.Label(demText));
  });
  panel.widgets().set(1, ui.Label(location));
// Edit: To be temporary, the "loading..." panel number has to be the same as the demText panel number (changed from 1 to 2).
  //panel.widgets().set(2, ui.Label("loading..."));
  Map.layers().set(1, ui.Map.Layer(click_point, {color: 'FF0000'}));
});

Geologi

  E so anche espatriare senza praticamente toccare strada asfaltata