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

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

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)

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

  );

Calcolo distanza tra 2 Aruco tags con OpenCV

Invece di calcolare la distanza tra la camera ed un tag in questo caso volevo vedere se riuscivo a calcolare la distanza tra due tag 

Raffronto la distanza reale tra due tag e la distanza calcolata da OpenCV. L'errore e' dell'ordine del cm

Nel video effettuato con un drone si vede come la distanza tra i due tag venga correttamente letta come costante nonostante il punto di riprese si muova

'''

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):

tvec0 = tvec

print(tvec0)

else:

tvec1 = tvec

print(tvec1)

distanza= "Dist=%4.0f"%(np.linalg.norm(tvec1-tvec0))

cv2.putText(frame, distanza,(50, 100),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()

Distanza e Pose Estimation da single aruco marker con opencv

Volevo vedere se potevo creare un distanziametro usando una webcam economica e gli aruco tag

Per le prove ho usato 

1) Realsense D435 ottico (1920x1080 2 MPx)
2) Logitech C920 (2560x1470 3.7 Mpx)
3) webcam cinese senza marca (3624 x 2448 8.8 Mpx)
4) DJI Tello video (1280x720 1Mpx)

Per la calibrazione ho utilizzato 5x5 e 9x6 su fogli A4

l'errore stimato di calibrazione e' risultato
Realsense : 0.0371
Logitech : 0.052
Cinese : 0.047
Tello video : 0.027

In generale maggiore e' la risoluzione della camera e maggiore e' la dimensione del tag piu' distante il tag viene riconosciuto. Con un tag da 25 cm 4x4 e la camera a maggiore risoluzione opencv riesce a vederlo fino a 30 m.  Con le altre camere non sono andato oltre i 20 m di distanza

Usando un telemetro laser ho calcolato la distanza del tag dalla camera e con Opencv ho calcolato la distanza. Come si vede dal grafico OpenCV riesce a risolvere la distanza .. da una prima stima c'e' un errore di 3-4 cm tra due posizioni vicine

prendendo solo i primi 10 metri la precisione aumenta

Ho provato anche a mettere a confronto gli Apriltag con gli Arucotag ma non ho notato sensibili vantaggi nell'usare AprilTags

Ho trovato invece miglioramenti nel raffinamento subpixel tra i parametri di OpenCV

elaborazione singole immagini

----------------------------------------------------------------

'''

Sample Command:-

python3 detect_aruco_images.py --image foto/103.png --type DICT_5X5_100 -a 0.1 -k ./calibration_matrix.npy -d ./distortion_coefficients.npy

'''

from turtle import width

import numpy as np

from utils import ARUCO_DICT, aruco_display

import argparse

import cv2

import sys

import math

ap = argparse.ArgumentParser()

ap.add_argument("-i", "--image", required=True, help="path to input image containing ArUCo tag")

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

#la dimensione del tag e' quella dello spigolo esterno del quadrato nero esterno, non i singoli quadrati interni

args = vars(ap.parse_args())

image = cv2.imread(args["image"])

print("Loading image "+ args["image"])

height,width,_ = image.shape

# don't resize for better perfomance

#h,w,_ = image.shape

#width=600

#height = int(width*(h/w))

#image = cv2.resize(image, (width, height), interpolation=cv2.INTER_CUBIC)

# verify that the supplied ArUCo tag exists and is supported by OpenCV

if ARUCO_DICT.get(args["type"], None) is None:

print(f"ArUCo tag type '{args['type']}' is not supported")

sys.exit(0)

# load the ArUCo dictionary, grab the ArUCo parameters, and detect

# the markers

#print("Detecting '{}' tags....".format(args["type"]))

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)

image = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)

arucoParams = cv2.aruco.DetectorParameters_create()

#--- 180 deg rotation matrix around the x axis

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

corners, ids, rejected = cv2.aruco.detectMarkers(image, arucoDict,parameters=arucoParams,cameraMatrix=k,distCoeff=d)

if len(corners) > 0:

#print(corners) #posizione degli angoli del marker

for i in range(0, len(ids)):

rvec, tvec, markerPoints = cv2.aruco.estimatePoseSingleMarkers(corners[i], float(args["aruco_dim"]), k,d)

#print(rvec)

dist = math.sqrt((tvec[0][0][0]*tvec[0][0][0])+(tvec[0][0][1]*tvec[0][0][1])+(tvec[0][0][2]*tvec[0][0][2]))

#dist = (dist*1.060672)+10.1674

print("Marker "+str(ids[i])+"="+str(dist))

-----------------------------------------------------------------

elaborazione video

-----------------------------------------------------------------

'''

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", required=True, 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('udp://127.0.0.1:5000',cv2.CAP_FFMPEG)

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.1)

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

for i in range(0, len(ids)):

rvec, tvec, markerPoints = cv2.aruco.estimatePoseSingleMarkers(corners[i], float(args["aruco_dim"]), k,d)

dist = math.sqrt((tvec[0][0][0]*tvec[0][0][0])+(tvec[0][0][1]*tvec[0][0][1])+(tvec[0][0][2]*tvec[0][0][2]))

distanza = "Distanza=%4.0f"%(dist)

cv2.putText(frame, distanza,(50, 50),cv2.FONT_HERSHEY_SIMPLEX, 1,(0, 0, 255), 2, cv2.LINE_4)

R_ct = np.matrix(cv2.Rodrigues(rvec)[0])

R_ct = R_ct.T

roll_marker, pitch_marker, yaw_marker = rotationMatrixToEulerAngles(R_flip*R_ct)

str_roll = "Roll=%4.0f"%(math.degrees(roll_marker))

str_pitch = "Pitch=%4.0f"%(math.degrees(pitch_marker))

str_yaw= "Yaw=%4.0f"%(math.degrees(yaw_marker))

cv2.putText(frame, str_roll,(50, 100),cv2.FONT_HERSHEY_SIMPLEX, 1,(0, 0, 255), 2, cv2.LINE_4)

cv2.putText(frame, str_pitch,(50, 125),cv2.FONT_HERSHEY_SIMPLEX, 1,(0, 0, 255), 2, cv2.LINE_4)

cv2.putText(frame, str_yaw,(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()

-----------------------------------------------------------------