Mascherare le nuvole in Sentinel 2

Nei dati Sentinel 2 e' contenuta anche in QI_DATA/MSK_CLDPRB_20m.jp2 una maschera di probabilita' che il pixel sia relativo ad una nuova (con valore di probabilita' crescente)

Per fare una maschera ho provato a settare i valori maggiori di 50 a zero ed i valori inferiori a 1 per poi moltiplicare i valori della banda

import numpy as np

import rasterio

from rasterio import plot

band4 = rasterio.open('b4.jp2', driver='JP2OpenJPEG')

mask = rasterio.open('mask.jp2', driver='JP2OpenJPEG')

band4_ = band4.read(1)

mask_ = mask.read(1)

mask_[mask_ <= 50] = 1

mask_[mask_ > 50] = 0

ris = np.zeros(band4.shape, dtype=rasterio.uint32)

ris =band4_*mask_

kwargs = band4.meta

kwargs.update(driver='GTiff',dtype=rasterio.uint16,count=1)

maskimage = rasterio.open('maschera0_1.tiff','w',**kwargs )

maskimage.write(mask_,1)

maskimage.close()

b4mask = rasterio.open('b4mask.tiff','w', **kwargs)

b4mask.write(ris,1)

b4mask.close()

Tagliare tutte le bande Sentinel in un colpo solo

 Puo' essere utile tagliare un granule di Sentinel sulla base di uno shape file

Per fare cio' si prepara lo shape (io ho usato lo stesso sistema di riferimento dell'immagine Sentinel... non so se rasterio funziona anche con sistema di riferimento misti) 

Tutte le bande sono nel folder immagini mentre il risultato sara' salvato nel folder rst

Nonostante sia Python e' uno script molto veloce

import os

import fiona

import rasterio

from rasterio.mask import mask

bandPath = './immagini/'

bandNames = os.listdir(bandPath)

aoiFile = fiona.open('./taglio.shp')

aoiGeom = [aoiFile[0]['geometry']]

for band in bandNames:

rasterPath = os.path.join(bandPath,band)

rasterBand = rasterio.open(rasterPath)

outImage, outTransform = mask(rasterBand, aoiGeom, crop=True)

outMeta = rasterBand.meta

outMeta.update({"driver": 'JP2OpenJPEG',

"height": outImage.shape[1],

"width": outImage.shape[2],

"transform": outTransform})

outPath = os.path.join('./rst',band)

outRaster = rasterio.open(outPath, "w", **outMeta)

outRaster.write(outImage)

outRaster.close()

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

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

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