sabato 12 novembre 2022

NDVI Sentinel 2 con RasterIO

il formato di output e' stato modificato in Geotiff perche' JP2 non permette di salvare in float32 


import numpy as np
import rasterio
from rasterio import plot

imagePath = './immagini/'
band4 = rasterio.open(imagePath+'B4.jp2', driver='JP2OpenJPEG')
band8 = rasterio.open(imagePath+'B8.jp2', driver='JP2OpenJPEG')

band4s = band4.read(1)
band8s = band8.read(1)

RED=band4s.astype(float)
NIR=band8s.astype(float)
np.seterr(divide = "ignore", invalid = "ignore")

ndvi = np.zeros(RED.shape, dtype=rasterio.float32)
ndvi = (NIR-RED) / (NIR+RED)


ndvi[ndvi > 1] = np.nan
ndvi[ndvi < -1] = np.nan

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

ndviimage = rasterio.open('./immagini/SentinelNDVI2.tiff','w',**kwargs )
ndviimage.write(ndvi,1)
ndviimage.close()


Truecolor Sentinel2 con RasterIO

Truecolor da dati Sentinel 2 usando RasterIO 


L'immagine e' stata corretta nel contrasto per renderla piu' leggibile


import rasterio
from rasterio import plot

imagePath = './immagini/'
band2 = rasterio.open(imagePath+'B2.jp2', driver='JP2OpenJPEG') #blue
band3 = rasterio.open(imagePath+'B3.jp2', driver='JP2OpenJPEG') #green
band4 = rasterio.open(imagePath+'B4.jp2', driver='JP2OpenJPEG') #red

band2s = band2.read(1)
band3s = band3.read(1)
band4s = band4.read(1)

bandg2=2.5*(band2s.astype(float))
bandg3=2.5*(band3s.astype(float))
bandg4=2.5*(band4s.astype(float))

trueColor = rasterio.open('./immagini/SentinelTrueColor2.tiff','w',driver='Gtiff',
width=band4.width, height=band4.height,
count=3,
crs=band4.crs,
transform=band4.transform,
dtype=band4.dtypes[0]
)
trueColor.write(bandg2,3) #blue
trueColor.write(bandg3,2) #green
trueColor.write(bandg4,1) #red
trueColor.close()

venerdì 11 novembre 2022

Sentinel 2 da SciHub

Nel precedente post avevo provato a scaricare i dati Sentinel da un bucket Google Cloud
In realta' ho trovato una procedura piu' lineare usando OpenHub di https://scihub.copernicus.eu/ 
E' necessario un account ma non e' legato a nessuna carta di credito ed il download e' libero




Per effettuare la selezione ed il download si usa lo script dhusget (Altre istruzioni a questo link )

una stringa di richiesta e' del tipo

./dhusget_0.3.8.sh -u username -p password -m Sentinel-2 -T S2MSI2A -t96 -o product -c 10.6,43.5:11.3,44 -n 1 -w 5 -O ./immagini

in cui viene effettuata la richeista delle immagini che contengano Firenze (vedi coordinate) acquisite nelle ultime 96 ore della missione Sentinel 2
Si possono specificare anche intervalli temporali tramite intervalli di date (in questo caso tra il 10/8/2022 ed il 12/8/2022

./dhusget_0.3.8.sh -u c1p81 -p xxxxxx@xxxxx -m Sentinel-2 -T S2MSI2A -S 2022-08-10T00:00:00.000Z -E 2022-08-12T00:00:00.000Z -n 1 -w 5 -c 10.6,43.5:11.3,44 -o product  -O ./immagini


Si possono scaricare sia il livello L1C che L2A .. di seguito si riporta il contenuto del file zip
I file zip vanno da circa 800 Mb ad oltre 1.1 Gb
Per selezionare solo un prodotto si puo' usare lo switch 
-T S2MSI1C per L21C
-T S2MSI2A per L2A

Ho notato che il sistema spesso tronca la connessione e non termina in modo corretto il download...se si abilita la funzione di checksum al termine non funziona mai..disabilitandolo invece si ottengono dei download corretti

LEVEL 1C
tutte le bande ognuna alla sua risoluzione nativa + truecolor a 10 m

LEVEL 2A 
risoluzione 10 m
AOT : Aerosol Optical Thickness
WVP : Water vapor 
TCI : Truecolor
Bande B2,B3,B4,B8

risoluzione 20 m
AOT : Aerosol Optical Thickness
WVP : Water vapor 
TCI : Truecolor
SCL : Scene classification map
Bande B1,B2,B3,B4,B5,B6,B7,B11,B12,B8A

risoluzione 60 m
AOT : Aerosol Optical Thickness
WVP : Water vapor 
TCI : Truecolor
SCL : Scene classification map
Bande B1,B2,B3,B4,B5,B6,B7,B9,B11,B12,B8A

con il medesimo sistema si possono ottenere i dati di Sentinel 1 SLC, OCN, GRDH e RAW

giovedì 10 novembre 2022

Sentinel 2 in Google Cloud

 Le immagini Sentinel2 sono rilasciate in modo gratuito ma i repository delle immagini le forniscono con una interfaccia utente interattiva mentre e' a pagamento l'accesso automatico via API

In modo abbastanza casuale mi sono imbattuto in questo bucket su Google Cloud

https://console.cloud.google.com/storage/browser/gcp-public-data-sentinel-2/tiles/33/U/UP;tab=objects?prefix=&forceOnObjectsSortingFiltering=false

che puo' essere utilizzato mediante questa libreria Python (per la selezione e download)

https://github.com/QuantuMobileSoftware/sentinel2tools

per l'installazione si usa

pip3 install 
git+https://github.com/QuantuMobileSoftware/sentinel2tools.git@ca232cb66106db6cac729defdab91aad9aecb15b.

Lo scripthe puo' essere utilizzato mediante questa libreria Python (per la selezione e download)

Il criterio di selezione geografica avviene tramite un geojson che puo' essere creato mediante il servizio on line geojson.io

from sentinel2download.downloader import Sentinel2Downloader
from sentinel2download.overlap import Sentinel2Overlap

if __name__ == '__main__':
verbose = True
aoi_path = "./firenze_geojson_utm.json"

overlap = Sentinel2Overlap(aoi_path, verbose=verbose)
tiles = overlap.overlap()

print(f"Overlapped tiles: {tiles}")

api_key = f"./api_key.json"

loader = Sentinel2Downloader(api_key, verbose=verbose)

product_type = 'L2A' # or L1C
start_date = "2022-09-01"
end_date = "2022-09-30"
output_dir = './utm'
cores = 3
BANDS = {'TCI', 'B04', }
CONSTRAINTS = {'NODATA_PIXEL_PERCENTAGE': 15.0, 'CLOUDY_PIXEL_PERCENTAGE': 10.0, }

loaded = loader.download(product_type,
tiles,
start_date=start_date,
end_date=end_date,
output_dir=output_dir,
cores=cores,
bands=BANDS,
constraints=CONSTRAINTS)

print(f"Load information")
for item in loaded:
print(item)

print("Execution ended")

Visto che i dati sono in un bucket di Gooogle Cloud e' richiesta anche 'autenticazione a GCloud. 

Per avere il file json di autenticazione si deve aprire la consolle di GCloud, creare un progetto e seguire le istruzioni al link seguente

https://cloud.google.com/api-keys/docs/get-started-api-keys

Le immagini vengono scaricate in formato JP2, una per banda, a 16 bit e georiferite in UTM (attenzione a questo dato). Oltre alle 12 bande sono presenti un truecolor (TCI) 




Viene scaricata tutta l'immagine che comprende la selezione geografica. Per effettuare il taglio si possono usare le GDAL. Le coordinate devono essere UTM

from osgeo import gdal

upper_left_x = 674882
upper_left_y = 4854678
lower_right_x = 688882
lower_right_y = 4844623
window = (upper_left_x,upper_left_y,lower_right_x,lower_right_y)

gdal.Translate('output_tci.jp2', 'tci.jp2', projWin = window)





Estrazione di piani da pointcloud

 Avevo gia' provato a fare una estrazione dei piani su un affioramento roccioso con PCL qualche tempo fa

Questa volta ho provato ad utilizzare la libreria Open3D con l'algoritmo Ransac per l'estrazione di piani multipli da una scansione fatta con il tablet Tango 



Per l'algoritmo e' sufficiente definire il numero massimo di piani attesi, la distanza di soglia ed il numero di iterazioni

I vari piani oltre ad essere visualizzati vengono anche estratti come ply

import sys
import numpy as np
import open3d as o3d
import matplotlib.pyplot as plt
input_path="./"
dataname="maiano.ply"
pcd = o3d.io.read_point_cloud(input_path+dataname)
o3d.visualization.draw_geometries([pcd])

#ransac piani multipli
segment_models={}
segments={}
max_plane_idx=3
rest=pcd
for i in range(max_plane_idx):
    colors = plt.get_cmap("tab20")(i)
    segment_models[i], inliers = rest.segment_plane(distance_threshold=0.1,ransac_n=3,num_iterations=5000)
    segments[i]=rest.select_by_index(inliers)
    segments[i].paint_uniform_color(list(colors[:3]))
    rest = rest.select_by_index(inliers, invert=True)
    print("pass",i,"/",max_plane_idx,"done.")
o3d.io.write_point_cloud("0.ply", segments[0])
o3d.visualization.draw_geometries([segments[0]] + [rest])
o3d.io.write_point_cloud("1.ply", segments[1])
o3d.visualization.draw_geometries([segments[1]] + [rest])
o3d.io.write_point_cloud("2.ply", segments[2])
o3d.visualization.draw_geometries([segments[2]] + [rest])

martedì 1 novembre 2022

Influenza dell'angolo di ripresa su Aruco Tags

 Usando gli Aruco Tags mi sono accorto che se vario la posizione di ripresa mantenendo i tag stabili gli angoli e le distanze fornite da RVec e TVec cambiano sensibilmente

Su github ho trovato questa libreria https://github.com/tentone/aruco in cui e' contenuta una tabella riassuntiva in cui si vede chiaramente l'effetto dell'angolo di ripresa sulla stima della distanza con errori anche sensibili per modeste (inferiori ai 10 gradi) variazioni dell'angolo di ripresa



in questo video si osserva che la distanza relativa tra i tag varia da un minimo di 31 ad un massimo di 40 cm con angoli tra 70 ed 80 al variare del punto di ripresa

sabato 29 ottobre 2022

StreamDesk DIY ESP32

 Tempo fa avevo gia' fatto un sistema a pulsanti per comandare Youtube...questa volta e' la versione definitiva. Pulsantiera collegata ad un ESP32 alimentato da 4 batterie AAA che emula la tastiera





Per l'alimentazione il pacco batterie e' collegata al VIN che e' un ingresso per tensione non regolata. Attenzione, non tutte le schede ESP32 hanno VIN ed al suo posto si trova un VOut a 5 V collegato direttamente alla USB

Per rendere i collegamenti semplici i pulsanti sono collegati direttamente all'ingresso digitale dell'ESP32 utilizzando l'impostazione PULLUP nel codice




#include <BleKeyboard.h>
BleKeyboard bleKeyboard;

#define BUTTON4 4 
#define BUTTON5 5 
#define BUTTON21 21 
#define BUTTON23 23 
int lState4 = HIGH; 
int cState4;   
int lState5 = HIGH; 
int cState5;   
int lState21 = HIGH; 
int cState21;   
int lState23 = HIGH; 
int cState23;   
void setup() {
  Serial.begin(9600);
  pinMode(BUTTON4, INPUT_PULLUP);
  pinMode(BUTTON5, INPUT_PULLUP);
  pinMode(BUTTON21, INPUT_PULLUP);
  pinMode(BUTTON23, INPUT_PULLUP);
  bleKeyboard.begin();
  
}
void loop() {
  if (bleKeyboard.isConnected()) {
  // read the state of the switch/button:
  cState4 = digitalRead(BUTTON4);
  cState5 = digitalRead(BUTTON5);
  cState21 = digitalRead(BUTTON21);
  cState23 = digitalRead(BUTTON23);
  if(lState4 == LOW && cState4 == HIGH)
      {
      Serial.println("DS4");
      bleKeyboard.print("K");
      }
      
  if(lState5 == LOW && cState5 == HIGH)
      {
      Serial.println("DS5");
      bleKeyboard.print("J");
      }
  if(lState21 == LOW && cState21 == HIGH)
      {
      Serial.println("DS21");
      bleKeyboard.print("F");
      }
  if(lState23 == LOW && cState23 == HIGH)
      {
      Serial.println("DS23");
      bleKeyboard.write(KEY_LEFT_ARROW);
      }
  
  lState4 = cState4;
  lState5 = cState5;
  lState21 = cState21;
  lState23 = cState23;
  delay(10);
  }
  else
  {
    Serial.println("In attesa di connessione");
    delay(5000);
  }
}

Distanza ed angoli relativi tra due tag Aruco

 La distanza reale centro centro tra i due tags e' di 35 cm



'''
Sample Command:-
python detect_aruco_video.py --type DICT_5X5_100 --camera True
python detect_aruco_video.py --type DICT_5X5_100 --camera False --video test_video.mp4 -a 25 -k ./calibration_matrix.npy -d ./distortion_coefficients.npy
'''

from turtle import delay
import numpy as np
from utils import ARUCO_DICT, aruco_display
import argparse
import time
import cv2
import sys
import math
import time

def isRotationMatrix(R):
Rt = np.transpose(R)
shouldBeIdentity = np.dot(Rt, R)
I = np.identity(3, dtype=R.dtype)
n = np.linalg.norm(I - shouldBeIdentity)
return n < 1e-6

def rotationMatrixToEulerAngles(R):
assert (isRotationMatrix(R))

sy = math.sqrt(R[0, 0] * R[0, 0] + R[1, 0] * R[1, 0])

singular = sy < 1e-6

if not singular:
x = math.atan2(R[2, 1], R[2, 2])
y = math.atan2(-R[2, 0], sy)
z = math.atan2(R[1, 0], R[0, 0])
else:
x = math.atan2(-R[1, 2], R[1, 1])
y = math.atan2(-R[2, 0], sy)
z = 0

return np.array([x, y, z])

R_flip = np.zeros((3,3), dtype=np.float32)
R_flip[0,0] = 1.0
R_flip[1,1] =-1.0
R_flip[2,2] =-1.0

ap = argparse.ArgumentParser()
ap.add_argument("-i", "--camera", help="Set to True if using webcam")
ap.add_argument("-v", "--video", help="Path to the video file")
ap.add_argument("-t", "--type", type=str, default="DICT_ARUCO_ORIGINAL", help="Type of ArUCo tag to detect")
ap.add_argument("-k", "--K_Matrix", required=True, help="Path to calibration matrix (numpy file)")
ap.add_argument("-d", "--D_Coeff", required=True, help="Path to distortion coefficients (numpy file)")
ap.add_argument("-a", "--aruco_dim", required=True, help="ArUco tag dimension")
args = vars(ap.parse_args())


if args["camera"].lower() == "true":
video = cv2.VideoCapture(0)
time.sleep(2.0)
else:
if args["video"] is None:
print("[Error] Video file location is not provided")
sys.exit(1)

video = cv2.VideoCapture(args["video"])

if ARUCO_DICT.get(args["type"], None) is None:
print(f"ArUCo tag type '{args['type']}' is not supported")
sys.exit(0)

arucoDict = cv2.aruco.Dictionary_get(ARUCO_DICT[args["type"]])
calibration_matrix_path = args["K_Matrix"]
distortion_coefficients_path = args["D_Coeff"]
k = np.load(calibration_matrix_path)
d = np.load(distortion_coefficients_path)
arucoParams = cv2.aruco.DetectorParameters_create()

while True:
ret, frame = video.read()
time.sleep(0.05)
if ret is False:
break


h, w, _ = frame.shape

width=1000
height = int(width*(h/w))
frame = cv2.resize(frame, (width, height), interpolation=cv2.INTER_CUBIC)
corners, ids, rejected = cv2.aruco.detectMarkers(frame, arucoDict, parameters=arucoParams)

detected_markers = aruco_display(corners, ids, rejected, frame)
if len(corners) > 0:
#print(corners) #posizione degli angoli del marker
#print(len(ids))
if (len(ids) == 2):
for i in range(0, len(ids)):
rvec, tvec, markerPoints = cv2.aruco.estimatePoseSingleMarkers(corners[i], float(args["aruco_dim"]), k,d)
if (i == 0):
print("Vettore 1")
tvec0 = tvec
R_ct = np.matrix(cv2.Rodrigues(rvec)[0])
R_ct = R_ct.T
roll, pitch, yaw = rotationMatrixToEulerAngles(R_flip*R_ct)
A = ([math.cos(roll),math.cos(pitch),math.cos(yaw)])
print(tvec0)
print(A)
if (i==1):
print("Vettore 2")
tvec1 = tvec
#B = ([rvec[0][0][0],rvec[0][0][1],rvec[0][0][2]])
R_ct = np.matrix(cv2.Rodrigues(rvec)[0])
R_ct = R_ct.T
roll, pitch, yaw = rotationMatrixToEulerAngles(R_flip*R_ct)
B = ([math.cos(roll),math.cos(pitch),math.cos(yaw)])
print(tvec1)
print(B)

# primo metodo per il calcolo della distanza
tvec0_x = tvec0[0][0][0]
tvec0_y = tvec0[0][0][1]
tvec0_z = tvec0[0][0][2]
tvec1_x = tvec1[0][0][0]
tvec1_y = tvec1[0][0][1]
tvec1_z = tvec1[0][0][2]
dist1 = math.sqrt(pow((tvec0_x-tvec1_x),2)+pow((tvec0_y-tvec1_y),2)+pow((tvec0_z-tvec1_z),2))
distanza1= "Dist=%4.0f"%(dist1)
#secondo metodo per il calcolo della distanza
#
distanza= "Dist=%4.0f"%(np.linalg.norm(tvec1-tvec0))
cv2.putText(frame, distanza1,(50, 100),cv2.FONT_HERSHEY_SIMPLEX, 1,(0, 0, 255), 2, cv2.LINE_4)
dot_product = np.dot(A,B,out=None)
normA = (np.linalg.norm(A))
normB = (np.linalg.norm(B))
cos_angolo = dot_product/(normA*normB)
angolo_rad = np.arccos(cos_angolo)
angolo_deg = np.rad2deg(angolo_rad)
if (angolo_deg > 90):
angolo_deg = 180 - angolo_deg
ang = "Ang=%4.1f"%(angolo_deg)
cv2.putText(frame,ang ,(50, 150),cv2.FONT_HERSHEY_SIMPLEX, 1,(0, 0, 255), 2, cv2.LINE_4)


cv2.imshow("Image", detected_markers)
key = cv2.waitKey(1) & 0xFF
if key == ord("q"):
break

cv2.destroyAllWindows()
video.release()

mercoledì 26 ottobre 2022

Calcolo di angolo tra vettori con NumPy

Calcolo di angolo tra due vettori con numpy secondo la formuale sottostante (prodotto scalare tra i due vettori diviso il prodotto delle norme dei vettori)

https://www.wikihow.it/Calcolare-l%27Angolo-tra-Due-Vettori
https://www.wikihow.it/Calcolare-l%27Angolo-tra-Due-Vettori


import numpy as np
A = np.array([3,3])
B = np.array([2,1])
dot_product = np.dot(A,B)
normA = (np.linalg.norm(A))
normB = (np.linalg.norm(B))
cos_angolo = dot_product/(normA*normB)
angolo_rad = np.arccos(cos_angolo)
print(np.rad2deg(angolo_rad))

martedì 25 ottobre 2022

Eclisse solare 25 Ottobre

 Ovviamente giornata nuvolosa che ha permesso pochi minuti di osservazione e mai in condizioni ottimali





venerdì 2 settembre 2022

Visualizzazione file GRIB meteo

 Per il download automatico e la visualizzazione di dati di previsione meteo ho trovato disponibile anche su Linux il software XyGrib

Permette di scaricare modelli di previsioni da diverse fonti di dati selezionando i campi di interesse su una finestra geografica definita dall'utente ed anche di creare delle animazioni (tasto razzo)

Anche solo per la comodita' vale la pena provarlo




giovedì 1 settembre 2022

Apriltag

 Dopo questa prova con Aruco Tags ho provato un progetto con AprilTag basato sulla libreria basato sulla libreria https://github.com/reinhartmeg/AprilTag-Distance-Calculator modificato per renderlo compatibile con OpenCV4

https://github.com/c1p81/distance_apriltag



mercoledì 24 agosto 2022

Coregistrazione immagini ottiche

 Questo post illustra come e' stato risolto il problema della GIF sottostante


Le immagini ottiche sono riprese da una camera PTZ che nel corso del giorno si sposta su diverse inquadrature...cio' comporta che non ritorna mai nelle impostazione perfette del giorno precedente rendendo impossibile l'analisi al computer ma anche la semplice interpretazione umana di eventuali modifiche

Per risolvere il problema sono state applicate tecniche che riguardano sia la correzione geometrica che dello spazio colore

Per prima cosa i vari fotogrammi sono stati coregistrati in modo software con il seguente script Python

///////////////////////////////////////////////////////

import cv2
import numpy as np
import sys, getopt

file2 = sys.argv[1:][0]
#print(file2)

#print("luca")
#print("/home/luca/ptz/150.217.73.10/share/transi"+file2)
# Read the images to be aligned
im1 = cv2.imread("/home/luca/ptz/150/share/transi/1.jpg");
im2 = cv2.imread("/home/luca/ptz/150/share/transi/"+file2);


# Convert images to grayscale
im1_gray = cv2.cvtColor(im1,cv2.COLOR_BGR2GRAY)
im2_gray = cv2.cvtColor(im2,cv2.COLOR_BGR2GRAY)

# Find size of image1
sz = im1.shape

# Define the motion model
warp_mode = cv2.MOTION_TRANSLATION

# Define 2x3 or 3x3 matrices and initialize the matrix to identity
if warp_mode == cv2.MOTION_HOMOGRAPHY :
warp_matrix = np.eye(3, 3, dtype=np.float32)
else :
warp_matrix = np.eye(2, 3, dtype=np.float32)

# Specify the number of iterations.
number_of_iterations = 5000;

# Specify the threshold of the increment
# in the correlation coefficient between two iterations
termination_eps = 1e-10;

# Define termination criteria
criteria = (cv2.TERM_CRITERIA_EPS | cv2.TERM_CRITERIA_COUNT, number_of_iterations, termination_eps)

# Run the ECC algorithm. The results are stored in warp_matrix.
(cc, warp_matrix) = cv2.findTransformECC (im1_gray,im2_gray,warp_matrix, warp_mode, criteria)

if warp_mode == cv2.MOTION_HOMOGRAPHY :
# Use warpPerspective for Homography
im2_aligned = cv2.warpPerspective (im2, warp_matrix, (sz[1],sz[0]), flags=cv2.INTER_LINEAR + cv2.WARP_INVERSE_MAP)
else :
# Use warpAffine for Translation, Euclidean and Affine
im2_aligned = cv2.warpAffine(im2, warp_matrix, (sz[1],sz[0]), flags=cv2.INTER_LINEAR + cv2.WARP_INVERSE_MAP);

# Show final results
#cv2.imshow("Image 1", im1)
#cv2.imshow("Image 2", im2)
#cv2.imshow("Aligned Image 2", im2_aligned)
cv2.imwrite("/home/luca/ptz/150/share/transi/all_"+file2,im2_aligned)
print("allinea")
print(file2)

///////////////////////////////////////////////////////

per la registrazione dell'istogramma dei colori e' stato usato

///////////////////////////////////////////////////////

from skimage import exposure
import matplotlib.pyplot as plt
import argparse
import cv2
# construct the argument parser and parse the arguments
ap = argparse.ArgumentParser()
ap.add_argument("-s", "--source", required=True,help="numero")
args = vars(ap.parse_args())

src = cv2.imread("/home/luca/ptz/150/share/transi/all_"+str(args["source"])+".jpg")
ref = cv2.imread("/home/luca/ptz/150/share/transi/all_1.jpg")
# determine if we are performing multichannel histogram matching
# and then perform histogram matching itself
print("[INFO] performing histogram matching...")
multi = True if src.shape[-1] > 1 else False
matched = exposure.match_histograms(src, ref, multichannel=multi)
cv2.imwrite("/home/luca/ptz/150/share/transi/orfeo"+str(args["source"])+".jpg",matched)

///////////////////////////////////////////////////////

per finire le immagini sono passate all'algoritmo di Change Detection di Orfeo ToolBox

otbcli_MultivariateAlterationDetector -in1 ~/ptz/150/share/transi/orfeo$j.jpg -in2 /home/luca/ptz/150.217.73.10/share/transi/orfeo$prec.jpg -out /home/luca/ptz/150/share/transi/change$j.jpg


ed il risultato viene elaborato da uno script in Python per esaltare il contrasto della mappa di CD

import cv2
import numpy as np
import sys, getopt

file2 = sys.argv[1:][0]


img = cv2.imread("/home/luca/ptz/150/transi/change"+file2+".jpg", cv2.IMREAD_COLOR)

# normalize float versions
#norm_img1 = cv2.normalize(img, None, alpha=0, beta=1, norm_type=cv2.NORM_MINMAX, dtype=cv2.CV_32F)
norm_img2 = cv2.normalize(img, None, alpha=0, beta=1.2, norm_type=cv2.NORM_MINMAX, dtype=cv2.CV_32F)

# scale to uint8
#norm_img1 = (255*norm_img1).astype(np.uint8)
norm_img2 = np.clip(norm_img2, 0, 1)
norm_img2 = (255*norm_img2).astype(np.uint8)

# write normalized output images
cv2.imwrite("/home/luca/ptz/150/transi/change_stretch"+file2+".jpg",norm_img2)
#cv2.imwrite("zelda1_bm20_cm20_normalize2.jpg",norm_img2)




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'})
  );

LLama3 Anita

A seguito di questo post ho provato a vedere ho provato a vedere cosa accadeva ad utilizzare un modello specifico per la lingua italiana in...