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

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

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

Previsioni temperatura con GFS NOAA ed Arpege

Ho provato a mettere a confronto le previsioni meteo dei servizi NOAA GFS ed Arpege visti nei precedenti post per una localita' italiana confrontando i dati con dati di verita' a terra dati da stazione meteo. Per rendere piu' agevole l'analisi il download ed il pretrattamento dei dati in CSV e' stato automatizzato  

Script per NOAA GFS

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

dove="/home/luca/gfs/"
rm $dove*.grib > /dev/null 2>&1
rm $dove*.csv > /dev/null 2>&1
rm /var/www/html/luca/gfs/completo.csv > /dev/null 2>&1
mv $dove*.txt ${dove}backup/ > /dev/null 2>&1

today=`date +%Y%m%d`
echo $today
# 11.8537, 46.4297
llon="11.75"
llat="46.25"
ulon="11.99"
ulat="46.49"
for i in $(seq -f "%03g" 0 120)
do
  echo $i
  wget -q -O ${dove}f${i}.grib "https://nomads.ncep.noaa.gov/cgi-bin/filter_gfs_0p25.pl?file=gfs.t00z.pgrb2.0p25.f$i&lev_2_m_above_ground=on&var_TMP=on&subregion=&leftlon=${llon}&rightlon=${ulon}&toplat=${ulat}&bottomlat=${llat}&dir=%2Fgfs.$today%2F00%2Fatmos"
done
for i in $(seq -f "%03g" 123 3 384) #123 126 129
do
  echo $i
  wget -q -O ${dove}f${i}.grib "https://nomads.ncep.noaa.gov/cgi-bin/filter_gfs_0p25.pl?file=gfs.t00z.pgrb2.0p25.f$i&lev_2_m_above_ground=on&var_TMP=on&subregion=&leftlon=${llon}&rightlon=${ulon}&toplat=${ulat}&bottomlat=${llat}&dir=%2Fgfs.$today%2F00%2Fatmos"
done

for i in $(seq -f "%03g" 0 120)
do
    echo "Conversione grib${i}"
  ${dove}wgrib2 ${dove}f${i}.grib -csv ${dove}f${i}.csv
done
for i in $(seq -f "%03g" 123 3 384)
do
  echo "Conversione grib${i}"
  ${dove}wgrib2 ${dove}f${i}.grib -csv ${dove}f${i}.csv
done
for i in $(seq -f "%03g" 0 120)
do
  cat ${dove}f${i}.csv >> ${dove}${today}_completo.txt
done
for i in $(seq -f "%03g" 123 3 384)
do
  cat ${dove}f${i}.csv >> ${dove}${today}_completo.txt
done
#aggiunge l'intestazione al csv
sed '1 s/^/data1,data2,variabile,unita,lon,lat,valore\n/' ${dove}${today}_completo.txt > ${dove}${today}_completo_1.txt
cp ${dove}${today}_completo_1.txt /var/www/html/luca/gfs/completo.csv

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

Script per Arpege (il servizio necessita di una API Key nel campo Token della URL ...ho modificato per non renderla utilizzabile la mia)

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

lon="11.8"

lat="46.4"

 

wget -O ${dove}00_12.grib2 "http://dcpc-nwp.meteo.fr/services/PS_GetCache_DCPCPreviNum?token=__5yLVTdr-sGeHoPitnFc7TZ6MhBcJxuxxxxxxxxxxxxxxx&model=ARPEGE&grid=0.1&package=SP1&time=00H12H&referencetime=${today}T00:00:00Z&format=grib2" 

wget -O ${dove}13_24.grib2 "http://dcpc-nwp.meteo.fr/services/PS_GetCache_DCPCPreviNum?token=__5yLVTdr-sGeHoPitnFc7TZ6MhBcxxxxxxxxxxxxxxxxxx&model=ARPEGE&grid=0.1&package=SP1&time=13H24H&referencetime=${today}T00:00:00Z&format=grib2" 

wget -O ${dove}25_36.grib2 "http://dcpc-nwp.meteo.fr/services/PS_GetCache_DCPCPreviNum?token=__5yLVTdr-sGeHoPitnFc7TZ6MhBcxxxxxxxxxxxxxxxxxx&model=ARPEGE&grid=0.1&package=SP1&time=25H36H&referencetime=${today}T00:00:00Z&format=grib2" 

wget -O ${dove}37_48.grib2 "http://dcpc-nwp.meteo.fr/services/PS_GetCache_DCPCPreviNum?token=__5yLVTdr-sGeHoPitnFc7TZ6MhBcJxuSsoZxxxxxxxxxxxxxxx&model=ARPEGE&grid=0.1&package=SP1&time=37H48H&referencetime=${today}T00:00:00Z&format=grib2" 

wget -O ${dove}49_60.grib2 "http://dcpc-nwp.meteo.fr/services/PS_GetCache_DCPCPreviNum?token=__5yLVTdr-sGeHoPitnFc7TZ6MhBcJxuSsoxxxxxxxxxxxxxxxxx&model=ARPEGE&grid=0.1&package=SP1&time=49H60H&referencetime=${today}T00:00:00Z&format=grib2" 

wget -O ${dove}61_72.grib2 "http://dcpc-nwp.meteo.fr/services/PS_GetCache_DCPCPreviNum?token=__5yLVTdr-sGeHoPitnFc7TZ6MhBcJxuxxxxxxxxxxxxxxxxxxxx&model=ARPEGE&grid=0.1&package=SP1&time=61H72H&referencetime=${today}T00:00:00Z&format=grib2" 

wget -O ${dove}73_84.grib2 "http://dcpc-nwp.meteo.fr/services/PS_GetCache_DCPCPreviNum?token=__5yLVTdr-sGeHoPitnFc7TZ6MhBxxxxxxxxxxxxxxxxxxxxxxx&model=ARPEGE&grid=0.1&package=SP1&time=73H84H&referencetime=${today}T00:00:00Z&format=grib2" 

wget -O ${dove}85_96.grib2 "http://dcpc-nwp.meteo.fr/services/PS_GetCache_DCPCPreviNum?token=__5yLVTdr-sGeHoPitnFc7TZ6MhBcJxxxxxxxxxxxxxxxxxxxxxx&model=ARPEGE&grid=0.1&package=SP1&time=85H96H&referencetime=${today}T00:00:00Z&format=grib2" 

wget -O ${dove}97_102.grib2 "http://dcpc-nwp.meteo.fr/services/PS_GetCache_DCPCPreviNum?token=__5yLVTdr-sGeHoPitnFc7TZ6MhBcJxxxxxxxxxxxxxxxxxxxxxxxx&model=ARPEGE&grid=0.1&package=SP1&time=97H102H&referencetime=${today}T00:00:00Z&format=grib2" 

${dove}wgrib2 -for 102:114 -s -lon ${lon} ${lat} ${dove}00_12.grib2 > ${dove}dati.csv

${dove}wgrib2 -for 97:108 -s -lon ${lon} ${lat} ${dove}13_24.grib2 >> ${dove}dati.csv

${dove}wgrib2 -for 97:108 -s -lon ${lon} ${lat} ${dove}25_36.grib2 >> ${dove}dati.csv

${dove}wgrib2 -for 97:108 -s -lon ${lon} ${lat} ${dove}37_48.grib2 >> ${dove}dati.csv

${dove}wgrib2 -for 97:108 -s -lon ${lon} ${lat} ${dove}49_60.grib2 >> ${dove}dati.csv

${dove}wgrib2 -for 97:108 -s -lon ${lon} ${lat} ${dove}61_72.grib2 >> ${dove}dati.csv

${dove}wgrib2 -for 97:108 -s -lon ${lon} ${lat} ${dove}73_84.grib2 >> ${dove}dati.csv

${dove}wgrib2 -for 97:108 -s -lon ${lon} ${lat} ${dove}85_96.grib2 >> ${dove}dati.csv

${dove}wgrib2 -for 97:108 -s -lon ${lon} ${lat} ${dove}97_102.grib2 >> ${dove}dati.csv

cut -d, -f1 --complement ${dove}dati.csv > ${dove}dati2.csv

cut -d, -f1 --complement ${dove}dati2.csv > ${dove}temp.csv

sed -i -e 's/val=//g' ${dove}temp.csv

go run ${dove}arpege.go > /var/www/html/luca/gfs/arpege.csv

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

Confronto tra GFS NOAA e stazione meteo a terra

Confronto tra Arpege e stazione meteo a terra

Come si vede in valore assoluto i dati sono differenti ma e' comunque possibile creare una retta di correlazione tra i dati di previsione da terra ed i dati misurati dalla stazione meteo. Considerando che il primo dato e' un dato di previsione il coefficiente di correlazione non e' pessimo

Evento Stromboli 12 agosto su ERA5

 Il 12 agosto 2022 e' avvenuto un inteso evento meteo a Stromboli ..volevo vedere se era possibile ricavare i dati di precipitazione anche dal servizio ERA5 di Copernicus per confrontarli con la verita' a terra data dal pluviometro INGV

per scaricare i dati ho utilizzato il seguente script Python (dati del solo 12 agosto per il solo parametro Total Precipitation)

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

import cdsapi

c = cdsapi.Client()

c.retrieve(

'reanalysis-era5-single-levels',

{

'product_type': 'reanalysis',

'variable': 'total_precipitation',

'year': '2022',

'month': '08',

'day': '12',

'time': [

'00:00', '01:00', '02:00',

'03:00', '04:00', '05:00',

'06:00', '07:00', '08:00',

'09:00', '10:00', '11:00',

'12:00', '13:00', '14:00',

'15:00', '16:00', '17:00',

'18:00', '19:00', '20:00',

'21:00', '22:00', '23:00',

],

'area': [

39, 15.25, 38.75,

15.5,

],

'format': 'netcdf',

},

'12.nc')

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

i dati sono stati acquisiti in formato NetCDF4 perche' il formato GRIB non sembra compatibile con Wgrib2 (forse si tratta di grib1)..per leggere i dati ho prima usato Qgis per verificare di avere centrato la finestra geografica di ricerca

e poi ho utilizzato il seguente script Python per estrarre i dati (non ho trovato niente di gia' pronto per una conversione da NetCDF4 a testo)

import netCDF4

precip_nc_file = '12.nc'

nc = netCDF4.Dataset(precip_nc_file, mode='r')

nc.variables.keys()

lat = nc.variables['latitude'][:]

lon = nc.variables['longitude'][:]

time_var = nc.variables['time']

dtime = netCDF4.num2date(time_var[:],time_var.units)

precip = nc.variables['tp'][:]

print("lat "+str(len(lat)))

print("lon "+str(len(lon)))

print("prec "+str(len(precip)))

print("time "+str(len(dtime)))

# coordinate Stromboli 38.79184084355282, 15.2152918018263

for s in range(len(dtime)):

for l in range(len(lat)):

for g in range(len(lon)):

print(dtime[s], end=''),

print(";", end='')

print(lat[l], end=''),

print(";", end='')

print(lon[g], end=''),

print(";", end='')

# il primo indice e' quello dell'ora

# il secondo indice e' la latitudine

# il terzo indice e' la longitudine

print(precip[s][l][g], end='\n\r')

questo e' il grafico generato con i dati di ERA5

questo per confronto sono i dati dal pluviometro di Stromboli gestito da INGV

(la scala dei tempi e' differente ed anche l'ora risulta differente a causa dei differenti fusi)

In conclusione il dato da satellite sottostima l'evento ma la forma generale del grafico e' raffrontabile

Debian e Ideapad 310 MIIX310

Mi sono comprato il tablet/convertibile Lenovo Ideapad 310 MIIX non perche' ne avessi bisogno ma per il costo estremamente ridotto (50 euro) e .... mi sono messo nei guai

Con Windows 10 la macchina e' completamente inutilizzabile ma me lo aspettavo....cio' che non aspettavo era di dover lottare con i driver per Linux

Una sintesi

con Ubuntu la scheda di rete Realtek8723BS funziona, non funziona audio e webcam (la webcam e' un modello molto particolare non USB ne' PCI ..e' una OV2680..comoda su Arduino per niente su PC per mancanza di supporto)

con Manjaro in versione live funziona praticamente tutto (tranne webcam)...una volta installato e' sparita la scheda di rete

Debian nemmeno a parlarne...il supporto alla Realtek 8723BS esisteva nel kernel fino a qualche anno fa ma e' stato rimosso. Ho provato ad installare Devuan per avere i driver gia' installati ...funziona l'audio ma niente wifi e webcam

alla fine pero' Debian e' molto rapida su una macchina dalle risorse cosi' limitate e mi sono adattato con una scheda Wifi USB

Un aspetto fastidioso e' che essendo un tablet il sistema operativo parte in modalita' portrait, per ruotare lo schermo e il touchscreen si possono utilizzare i comani

xrandr -o right   
xinput set-prop 'FTSC1000:00 2808:1015' 'Coordinate Transformation Matrix' 0 1 0 -1 0 1 0 0 1
 

 

 sdfsfd

Distanza Realsense

Misura la distanza al centro dell'immagine

Per rendere la misura piu' stabile e' stato utilizzato un ring buffer

## License: Apache 2.0. See LICENSE file in root directory.

## Copyright(c) 2015-2017 Intel Corporation. All Rights Reserved.

###############################################

## Open CV and Numpy integration ##

###############################################

import pyrealsense2 as rs

import numpy as np

from numpy_ringbuffer import RingBuffer

import cv2

RingBufferSize = 50

r = RingBuffer(capacity=RingBufferSize, dtype=np.float16)

# Configure depth and color streams

pipeline = rs.pipeline()

config = rs.config()

# Get device product line for setting a supporting resolution

pipeline_wrapper = rs.pipeline_wrapper(pipeline)

pipeline_profile = config.resolve(pipeline_wrapper)

device = pipeline_profile.get_device()

device_product_line = str(device.get_info(rs.camera_info.product_line))

found_rgb = False

for s in device.sensors:

if s.get_info(rs.camera_info.name) == 'RGB Camera':

found_rgb = True

break

if not found_rgb:

print("The demo requires Depth camera with Color sensor")

exit(0)

config.enable_stream(rs.stream.depth, 640, 480, rs.format.z16, 30)

if device_product_line == 'L500':

config.enable_stream(rs.stream.color, 960, 540, rs.format.bgr8, 30)

else:

config.enable_stream(rs.stream.color, 640, 480, rs.format.bgr8, 30)

# Start streaming

pipeline.start(config)

try:

while True:

# Wait for a coherent pair of frames: depth and color

frames = pipeline.wait_for_frames()

depth_frame = frames.get_depth_frame()

color_frame = frames.get_color_frame()

if not depth_frame or not color_frame:

continue

# Convert images to numpy arrays

depth_image = np.asanyarray(depth_frame.get_data())

color_image = np.asanyarray(color_frame.get_data())

#print(depth_image[320,240])

r.append(depth_image[320,240])

print(np.mean(r))

# Apply colormap on depth image (image must be converted to 8-bit per pixel first)

depth_colormap = cv2.applyColorMap(cv2.convertScaleAbs(depth_image, alpha=0.03), cv2.COLORMAP_JET)

depth_colormap_dim = depth_colormap.shape

color_colormap_dim = color_image.shape

# If depth and color resolutions are different, resize color image to match depth image for display

if depth_colormap_dim != color_colormap_dim:

resized_color_image = cv2.resize(color_image, dsize=(depth_colormap_dim[1], depth_colormap_dim[0]), interpolation=cv2.INTER_AREA)

images = np.hstack((resized_color_image, depth_colormap))

else:

images = np.hstack((color_image, depth_colormap))

# Show images

cv2.namedWindow('RealSense', cv2.WINDOW_AUTOSIZE)

startpoint=(310,240)

endpoint=(330,240)

color=(255,0,0)

thickness=3

images = cv2.line(images,startpoint,endpoint,color,thickness)

startpoint=(320,230)

endpoint=(320,250)

images = cv2.line(images,startpoint,endpoint,color,thickness)

cv2.imshow('RealSense', images)

cv2.waitKey(1)

finally:

# Stop streaming

pipeline.stop()