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

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


martedì 18 aprile 2023

Cuda Mandelbrot

 


 

 

#include <stdio.h>
//questa libreria crea png senza altre dipendenze
#include "libattopng.h"

#define W  5000 //width image
#define H  5000 //height image
#define TX 32 // number of threads per block along x-axis
#define TY 32 // number of threads per block along y-axis

__global__
void Kernel(int *d_out, int w, int h, int itera)
{
    //nel kernel viene lanciato un thread per ogni pixel dell'immagine
    //per calcolare la riga e la colonna si usa la formula sottostante
    int c = blockIdx.x*blockDim.x + threadIdx.x; // colonna
    int r = blockIdx.y*blockDim.y + threadIdx.y; // riga
    int i = r*w + c; // posizione del pixel nell'array lineare
    if ((c >= w) || (r >= h)) return;

    // inizia il ca
    const float x_min = -0.727;
    const float x_max = -0.730;
    const float y_min = -0.247;
    const float y_max = -0.250;

    float dx,dy;

    dx = (x_max-x_min)/w;
    dy = (y_max-y_min)/h;


    float k = x_min+(r*dx);
    float j = y_min+(c*dy);
    d_out[i] = 0;

    double x,y,x_new, y_new = 0.0;
    for (int f=0; f<itera;f++)
            {
                    x_new = (x*x) - (y*y) + k;
                    y_new = (2*x*y) +j;
                    if (((x_new*x_new)+(y_new*y_new))>4)
                    {
                        d_out[i] = f;
                        return;
                    }
                    x = x_new;
                    y = y_new;
            }
 }

int main()
{
 // alloca e mette a zero (malloc alloca e basta) un array di int di dimensione u
 // uguale all'immagine desiderata sulla ram della CPU
 int *out = (int*)calloc(W*H, sizeof(int));

 // alloca un array nella GPU delle stesse dimensioni
 // la GPU viene definita come device mentre la CPU e' l'host
 int *d_out; // pointer for device array
 cudaMalloc(&d_out, W*H*sizeof(int));

 // definisce i threads
 const dim3 blockSize(TX, TY);
 const int bx = (W + TX - 1)/TX;
 const int by = (W + TY - 1)/TY;
 const dim3 gridSize = dim3(bx, by);
 //lancia il kernel sulla GPU passando come parametro l'array sulla GPU
 Kernel<<<gridSize, blockSize>>>(d_out, W, H,4096);

 //trasferisce il contenuto dell'array dalla GPU alla RAM della CPU
 cudaMemcpy(out, d_out, W*H*sizeof(int), cudaMemcpyDeviceToHost);

 //crea una PNG partendo dall'array
 libattopng_t* png = libattopng_new(W, H, PNG_GRAYSCALE);
 
 int x, y;
 for (y = 0; y < W; y++) {
   for (x = 0; x < H; x++) {
     libattopng_set_pixel(png, y, x, out[x+y*H]%255);
   }
 }
 libattopng_save(png, "mandelbrot.png");
 libattopng_destroy(png);


 cudaFree(d_out);
 free(out);
 return 0;
 }
 

 

mercoledì 12 aprile 2023

Trasformazione affine con OpenCV

 Alla fine degli anni 80 se non erro alla Domenica Sportiva c'era un  primordiale sistema di computer grafica che consentiva di valutare il fuorigioco in una partita di calcio da immagini televisive....mi ero sempre chiesto come era possivile....dopo tanto tempo ho visto che e' possibile tramite la trasformazione affine 

 Partiamo da una immagine generica dello stadio di Firenze


 Le dimensione del campo dell'Artemio Franchi sono 105x68 m mentre nell'immagine sottostante sono riportate le dimensioni ufficiali del campo da calcio


 La trasformazione affine serve a trasformare le coordinate immagine in un altro sistema di coordinate (in questo caso coordinate centrimetriche reali con origine degli assi nella bandierina del calcio d'angolo in basso a sinistra)
OpenCV richiede tre punti ...per esempio la bandierina del calcio d'angolo in alto a sinistra ha coordinate pixel di 362x208 mentre coordinate reali 10500,6800 cm. I punti sono stati selezionati sull'incrocio di alcune linee dell'area di rigore
Nell'immagine successiva il risultato dell'algoritmo...non perfetto perche' le linee dovrebbero essere ortogonali ma considerando le incertezze direi che non e'male

 
import matplotlib.pyplot as plt
import numpy as np
import cv2 as cv

img = cv.imread('fiorentina.jpg')
pts1 = np.float32([[193,132],[362,208],[542,214]])
pts2 = np.float32([[10500,6800],[9950,4315],[10500,3765]])

M = cv.getAffineTransform(pts1,pts2)
dst = cv.warpAffine(img,M,(10500,6800))

plt.imshow(dst)
plt.show()
plt.savefig("dest.png")
plt.close(dst)


 

mercoledì 29 marzo 2023

PostgresSql partitioning

Una delle funzioni meno comuni di Postgres e' quella di effettuara il partiniong di una tabella per velocizzare il tempo di accesso. In sintesi si puo' suddividere una tabella di generose dimensione sulla base di un criterio (data, indice numerico intero) in sotto tabelle senza modificare la query di interrogazione della tabella madre. Non tutte le strutture dati sono idonee al partizionamento (per esempio quando ci sono solo campi stringhe)

In pratica si puo' creare un partitioning inserendo il campo che si vuole utilizzare per il partitioning

CREATE TABLE public.partizione (

    code character varying(12),
    vel real,
    v_stdev real,
  ....
    objectid integer NOT NULL,
    the_geom public.geometry(Point,32633)
) PARTITION BY RANGE(objectid);

 si creano quindi le sottotabelle con gruppi di 20000....attenzione che se il valore del campo di partizionamento eccede il valore massimo la riga non verra' inserita...per esempio inserire un campo con objectid=290000 generera' errore con la regola sottostante

CREATE TABLE partizione_0 PARTITION OF partizione FOR VALUES FROM (1) TO (20000);
CREATE TABLE partizione_1 PARTITION OF partizione FOR VALUES FROM (20001) TO (40000);
CREATE TABLE partizione_2 PARTITION OF partizione FOR VALUES FROM (40001) TO (60000);
CREATE TABLE partizione_3 PARTITION OF partizione FOR VALUES FROM (60001) TO (80000);
CREATE TABLE partizione_4 PARTITION OF partizione FOR VALUES FROM (80001) TO (100000);
CREATE TABLE partizione_5 PARTITION OF partizione FOR VALUES FROM (100001) TO (120000);
CREATE TABLE partizione_6 PARTITION OF partizione FOR VALUES FROM (120001) TO (140000);
CREATE TABLE partizione_7 PARTITION OF partizione FOR VALUES FROM (140001) TO (160000);
CREATE TABLE partizione_8 PARTITION OF partizione FOR VALUES FROM (160001) TO (180000);
CREATE TABLE partizione_9 PARTITION OF partizione FOR VALUES FROM (180001) TO (200000);
CREATE TABLE partizione_10 PARTITION OF partizione FOR VALUES FROM (200001) TO (220000);
CREATE TABLE partizione_11 PARTITION OF partizione FOR VALUES FROM (220001) TO (240000);
CREATE TABLE partizione_12 PARTITION OF partizione FOR VALUES FROM (240001) TO (260000);
CREATE TABLE partizione_13 PARTITION OF partizione FOR VALUES FROM (260001) TO (280000);

Per interrogare la tabella si usa SELECT..FROM TABLE partizione...sara' poi il motore del DB ad andare a trovare la sottotabella corrispondente

Ho cercato di misurare l'aumento delle prestazione del database ma ho avuto risultati variabili anche a causa di non riuscire a controllare la cache...indicativamente il miglior risultato e' stata un diminuzione dei tempi del 30%

giovedì 16 marzo 2023

Phishing via SMS 2023

 Una variante interessante rispetto al precedente esempio di phishing via SMS visto in precedenza 

Questa volta e' stato piu' insidioso perche' viene indicata la consegna di un pacco ad un centro di ritiro...plausibile...chi e' che non riceve pacchi oggi


Andando al link viene segnalata la necessita' di pagare 2 euro di tasse doganali (plausibile...mi e' capitato di ricevere pacchi dall'estero con la necessita' di pagare le tasse)




Mi sono insospettito quando mi veniva richiesto di reinserire i miei dati personali (del resto il corriere li doveva gia' conoscere....ho guardato il sito awardau.live..e sono passato ad aprire il link sul computer per vedere i sorgenti della pagina. La cosa divertente e' andando a www.awardau.live veniva proposto una spedizione di DHL mentre il link dell'SMS indicava una spedizione GLS


Capito l'ingannno ho provato a vedere come funziona la truffa....dopo aver cliccato sul sito awardau.live si veniva instradati su Pagamento Sicuro (ovviamente finto) in cui vengono richieste tutte le informazioni sulla carta di credito (ovviamente non era presente il pagamento via Paypal) compreso il CIN.

Nel codice HTML oltre a riferimenti a JQuery era presente il link a https://track.fwdtrck.com/click
(il link e' disattivo al momento in cui scrivo) ed a pro.ip-api.con con una key scaduta





Pulsante Dyson V10

 Pagare oltre 400 euro per avere un aspirapolvere e dopo tre anni essere costretti a smontarlo per la rottura di pezzo di plastica (il pulsante) da pochi centesimi di euro non fa certo piacere....fra le altre cose per la riparazione si deve smontare praticamente tutto il dispositivo (e si deve almeno una chiave torx con prolunga)

Il pezzo di ricambio corregge la debolezza di progettazione iniziale





mercoledì 15 marzo 2023

Local Packages in Golang

 Per poter utilizzare packages locali (librerie di funzioni create dall'utente esterno racchiuse in un package esterno al main) la soluzione e' inserire i le librerie esterne nel subfolder /vendor



nel file main.go l'importazione avviene semplicemente inserendo in import il nome del file della libreria e richiamando le funzioni della libreria come metodi 


import (
    "fmt"
    "scrivi"
)

func main() {
    fmt.Println("luca")
    scrivi.Debug(true)
    scrivi.Log("log")
}


Debugger integrato ESP32S3

Aggiornamento In realta' il Jtag USB funziona anche sui moduli cinesi Il problema risiede  nell'ID USB della porta Jtag. Nel modulo...