RGB change detection con RasterIO

 Un semplice metodo per avere il change detection da una immagine RGB e' la formula seguente

 

 che non e' altro che la distanza euclidea tra i due punti rappresentativi delle immagini nello spazio delle bande (quindi puo' essere scalato anche su n-bande)

 

 utilizzando la libreria RasterIO si puo' scrivere il codice seguente

 

import numpy as np

import rasterio

from rasterio.plot import show_hist,show

immagine1 = "cd1.tif"

immagine2 = "cd2.tif"

im1 = rasterio.open(immagine1)

im2 = rasterio.open(immagine2)

im1_b2 = im1.read(1).astype(float)

im1_b3 = im1.read(2).astype(float)

im1_b4 = im1.read(3).astype(float)

im2_b2 = im2.read(1).astype(float)

im2_b3 = im2.read(2).astype(float)

im2_b4 = im2.read(3).astype(float)

np.seterr(divide = "ignore", invalid = "ignore")

change_det = np.zeros(im1_b2.shape, dtype=rasterio.float32)

change_det = np.sqrt(np.power(im1_b2-im2_b2,2)+np.power(im1_b3-im2_b3,2)+np.power(im1_b4-im2_b4,2))

#print("Change Det Max " + str(np.max(change_det)))

#print("Change Det Min " + str(np.min(change_det)))

#show(change_det)

#show_hist(change_det)

kwargs = im1.meta

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

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

cd_image.write(change_det,1)

cd_image.close()

 

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

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