venerdì 3 marzo 2023

Telit GE864-PY

Mi e' stato regalato una volta dismesso questo dispositivo nato per leggere una seriale ed inviare i dati via rete cellulare 


Pensavo di poter recuperare qualche componente ma tutte le funzioni sono integrate nell'unico chip



Seagate ST-251

 Il mio piccolo museo dell'informatica ha avuto una nuova acquisizione, un hard disk Seagate ST-251 del 1990

Si tratta di un HD da 42 Mb formattato che ha connettori MTM e non IDE. Non ho il controller adatto per cui non ho la piu' pallida idea se sia funzionante (ma guardando su Internet sembra che questi dispositivi fossero particolarmente delicati e pochi sono sopravvissuti funzionanti)







mercoledì 1 marzo 2023

PIP Python error: externally-managed-environment

 Da qualche tempo sulla mia Debian testing non riesco ad installare le librerie Python tramite PIP ma sono tramite APT a causa dell'errore

error: externally-managed-environment

× This environment is externally managed
╰─> To install Python packages system-wide, try apt install python3-xyz, where xyz is the package you are trying to install.


l'unico modo per risolvere e' quello di usare un venv od usare PyCharm che in automatico crea un venv per ogni progetto

Change detection con Open3D su nuvole di punti

Anche la libreria Open3D ha una funzione che permette di calcolare la distanza tra due nuvole di punti anche se non e' esplicitato l'algoritmo (inoltre per i punti in cui non c'e' sovrapposizione immediata il codice del precedente post riporta NaN mentre qui un valore di distanza viene sempre proposto. Inoltre la matrice delle distanze contiene solo valori positivi



I risultati mi sembrano migliori quelli dell'algoritmo del precedente post

import numpy as np
import laspy as lp
import open3d as o3d

p2019= lp.read("2019.las")

points_2019 = np.vstack((p2019.x, p2019.y, p2019.z)).transpose()
#colors_2019 = np.vstack((p2019.red, p2019.green,p2019.blue)).transpose()

pcd2019 = o3d.geometry.PointCloud()
pcd2019.points = o3d.utility.Vector3dVector(points_2019)
#pcd2019.colors = o3d.utility.Vector3dVector(colors_2019 / 65535)


p2022 = lp.read("2022.las")

points_2022 = np.vstack((p2022.x, p2022.y, p2022.z)).transpose()
#colors_2022 = np.vstack((p2022.red, p2022.green,p2022.blue)).transpose()

pcd2022 = o3d.geometry.PointCloud()
pcd2022.points = o3d.utility.Vector3dVector(points_2022)
distance=pcd2019.compute_point_cloud_distance(pcd2022)

rows, columns = points_2019.shape

reds = np.empty(rows)
blues = np.empty(rows)
greens = np.empty(rows)


for i in range(0,rows):
#print(t)
if (distance[i]< -1.5):
Red = 255
Green = 0
Blue = 0

if ((distance[i]>=-1.5) and (distance[i]<=-0.75)):
Red = 255
Green = 255
Blue = 0

if ((distance[i]>=-0.75) and (distance[i]<=0.75)):
Red = 0
Green = 255
Blue = 0
if ((distance[i]>=0.75) and (distance[i]<=1.5)):
Red = 0
Green = 255
Blue = 255

if (distance[i]> 1.5):
Red = 0
Green = 0
Blue = 255
if (distance[i]> 10):
Red = 0
Green = 0
Blue = 0

if (distance[i]< 10):
Red = 0
Green = 0
Blue = 0

reds[i]= Red
greens[i]= Green
blues[i]=Blue

header = lp.LasHeader(point_format=3, version="1.2")
las = lp.LasData(header)

las.x = points_2019[:,0]
las.y = points_2019[:,1]
las.z = points_2019[:,2]
las.red = reds[:]
las.blue = blues[:]
las.green = greens[:]
las.write("open3d_result.las")






lunedì 27 febbraio 2023

Change detection con M3C2 su nuvole di punti

L'algoritmo M3C2 viene utlizzato da Cloudcompare per effettuare il change detection ...volevo trovare una soluzione per rendere automatico il processing e la pubblicazione

Si parte da due LAS e tramite la libreria py4dgeo si ottiene un file las con i punti colorati secondo una scala colore basata sulla distanza tra i due LAS si origine



import numpy as np
import py4dgeo
import laspy
import sys
import math

py4dgeo.set_num_threads(4) # numero di threads attivi
epoch1, epoch2 = py4dgeo.read_from_las("2019.las", "2022.las")
corepoints = epoch1.cloud[::30]

m3c2 = py4dgeo.M3C2(epochs=(epoch1, epoch2),corepoints=corepoints,cyl_radii=(2.0,),normal_radii=(0.5, 1.0, 2.0),)
distances, uncertainties = m3c2.run()

dist_max = np.nanmax(distances)
dist_min = np.nanmin(distances)

reds = np.empty(corepoints.shape)
blues = np.empty(corepoints.shape)
greens = np.empty(corepoints.shape)

fattore =(dist_max-dist_min)/(dist_max+dist_min)

Red =0
Blue = 0
Green = 0
for i in range (distances.shape[0]):
if math.isnan(distances[i]):
t = 0
else:
#print(distances[i])
t = int(((distances[i]-dist_min)/(dist_max-dist_min)) * 16581375)
#print(t)

if math.isnan(t):
Blue = 0
Green = 0
Red = 0
else:
if (distances[i]< -1.5):
Red = 255
Green = 0
Blue = 0

if ((distances[i]>=-1.5) and (distances[i]<=-0.75)):
Red = 255
Green = 255
Blue = 0

if ((distances[i]>=-0.75) and (distances[i]<=0.75)):
Red = 0
Green = 255
Blue = 0
if ((distances[i]>=0.75) and (distances[i]<=1.5)):
Red = 0
Green = 255
Blue = 255

if (distances[i]> 1.5):
Red = 0
Green = 0
Blue = 255
#RGBint = t
#Red = RGBint & 255
#Blue = (RGBint >> 8) & 255
#Green = (RGBint >> 16) & 255'''



reds[i]= Red
greens[i]= Green
blues[i]=Blue

header = laspy.LasHeader(point_format=3, version="1.2")
header.add_extra_dim(laspy.ExtraBytesParams(name="distance", type=np.float32))

xmin = np.floor(np.min(corepoints[:,0]))
ymin = np.floor(np.min(corepoints[:,1]))
zmin = np.floor(np.min(corepoints[:,2]))

header.offsets = [xmin,ymin,zmin]
#header.scales = np.array([0.1, 0.1, 0.1])
las = laspy.LasData(header)

las.x = corepoints[:,0]
las.y = corepoints[:,1]
las.z = corepoints[:,2]
las.red = reds[:,0]
las.blue = blues[:,0]
las.green = greens[:,0]
las.distance = np.copy(distances)
las.write("m3c2_result.las")


#./PotreeConverter m3c2_result.las -o /var/www/html/m3c2/ --generate-page index


Per la pubblicazione su web il file viene convertito e vestito tramite PoTree 

 
  • Anders, K., Winiwarter, L., Lindenbergh, R., Williams, J. G., Vos, S. E., & Höfle, B. (2020). 4D objects-by-change: Spatiotemporal segmentation of geomorphic surface change from LiDAR time series. ISPRS Journal of Photogrammetry and Remote Sensing, 159, pp. 352-363. doi: 10.1016/j.isprsjprs.2019.11.025.

  • Anders, K., Winiwarter, L., Mara, H., Lindenbergh, R., Vos, S. E., & Höfle, B. (2021). Fully automatic spatiotemporal segmentation of 3D LiDAR time series for the extraction of natural surface changes. ISPRS Journal of Photogrammetry and Remote Sensing, 173, pp. 297-308. doi: 10.1016/j.isprsjprs.2021.01.015.

  • Truong, C., Oudre, L., Vayatis, N. (2018): ruptures: Change point detection in python. arXiv preprint: abs/1801.00826.

giovedì 23 febbraio 2023

Composit Sentinel 2

Per creare un composit da Snap con dati Sentinel 2 si deve fare prima un subset (Raster/Subset) delle tre bande richieste. Se la risoluzione delle bande e' differente si deve successivamente rendere la grandezza dei pixel omogenea con Geometric/Resampling

    431

A questo punto Open RGB Image Window....in Colour Management si puo' impostare lo stretch


venerdì 17 febbraio 2023

Primi passi con R PCA

Si parte da un file CSV in cui le colonne sono le variabili (in questo caso dati analitici) 

La prima riga e' di intestazione, separatore punto e virgola, punto decimake

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

Punto;As;Co;Hg
Pz3A;26.4;13.4;0.206
Pz3B;41.9;17.7;0.281
Pz4;8.2;6.05;1.13
Pz5A;21.5;18.6;0.281
Pz5B;20.1;18.2;0.235
Pz6;5.8;20.4;0.138
Pz7A;6.8;23.2;0.082
Pz7B;7.4;24;0.090
Pz8A;16.8;21.8;0.325
Pz8B;12.4;22.2;0.279
Pz9A;46.9;30.2;0.71
Pz9B;20.0;21.6;0.319
Pz10A;29.3;9.4;1.03
Pz10B;18.5;56.1;0.45
Pz11A;63;12.7;0.210
Pz12;32.1;3.53;0.0298
Pz13A;18.2;31.1;0.75
Pz13B;14.8;29.1;0.200
Pz14A;9.5;28.6;0.043
Pz14B;9.2;32.6;0.045

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

Per elaborare i dati con PCA vengono caricate le librerie FactoMineR e factoextra

Per leggere il file CSV si deve specificare il separatore decimale

Si deve anche indicare che la prima colonna e' un dato qualitativo (da escludere dal calcolo), le dimensioni della PCA pari a 3,




 

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

install.packages("FactoMineR")
install.packages("factoextra")

library("FactoMineR")
library("factoextra")
data_res <- read.csv2("c:/Users/l.innocenti/Desktop/arsenico/arsenico.csv",dec=".")
data_res.pca <- PCA(data,scale.unit = TRUE,ncp=3,graph = TRUE,quali.sup=1)
data_res.pca
eig.val <- get_eigenvalue(data_res.pca)
fviz_eig(data_res.pca,addlabels = TRUE)
fviz_pca_var(data_res.pca)
var <-get_pca_var(data_res.pca)
library(corrplot)
corrplot(var$cos2,is.corr = FALSE)
fviz_pca_ind(data_res.pca)

Red Edge con Sentinel 2

Volevo migliorare un po' quanto provato qui piu' che altro per avere una migliore risoluzione spaziale. Ho provato con Sentinel 2 (...