Introduzione Data Wrangling di Immagini Geospaziali 🛰️ 🌍¶
In questo tutorial vedremo come poter accedere ai dati del satellite Sentinel-2. Vedremo alcune librerie comuni utilizzate nel Geospatial Data Science ed il loro utilizzo per accere ai dati utili al progetto Copernicus.
In particolare vedremo:
- Come generare un footprint (ossia come definire l'area di studio)
- Come interrogare il database di delle immagini delle missioni Sentinel (ossia il Copernicus Open Access Hub)
- Come scaricare ed estrarre i file interessanti
- Come leggere, scrivere ed infine manipolare le immagini riprese.
Tutorial by Francesco Pelosin
Installazione prerequisiti¶
In [0]:
!apt install ffmpeg
!pip install geopandas
!pip install folium
!pip install sentinelsat
!pip install rasterio
!pip install opencv-python
!pip install imageio
In [0]:
# Librerie specifiche per lavorare con dati Geospaziali
import geopandas as gpd
import folium
from sentinelsat import SentinelAPI, read_geojson, geojson_to_wkt
# Librerie specifiche per lavorare con dati immagine
import rasterio as rio
import rasterio.mask
from rasterio.plot import show
import matplotlib.pyplot as plt
import imageio
import cv2
# Libreria calcolo vettoriale
import numpy as np
# Librerie di aiuto per l'interfaccia OS
import zipfile
import os
import base64
from IPython.display import HTML
In [0]:
# Preparazione delle cartelle necessarie
if not os.path.exists("./data/"):
os.makedirs("./data/")
if not os.path.exists("./video/"):
os.makedirs("./video/")
if not os.path.exists("./imgs/"):
os.makedirs("./imgs/")
if not os.path.exists("./geoj/"):
os.makedirs("./geoj/")
def rim_info(rim):
print(f"Number of bands: {rim.count}")
print(f"Dset bounds: {rim.bounds}")
print(f"{(rim.bounds.top - rim.bounds.bottom)/ 1000} km² high")
print(f"{(rim.bounds.right - rim.bounds.left)/ 1000} km² wide")
def playvideo(filename):
video = open(filename, 'r+b').read()
encoded = base64.b64encode(video)
return HTML(data='''<video alt="test" controls><source src="data:video/mp4;base64,{0}" type="video/mp4"/></video>'''.format(encoded.decode('ascii')))
def get_bands_files(R_dir):
bands_dict = {}
for x in os.listdir(R_dir):
if ('_B08_' in x):
bands_dict['B08'] = R_dir+x
elif ('_B12_' in x):
bands_dict['B12'] = R_dir+x
elif ('_B11_' in x):
bands_dict['B11'] = R_dir+x
elif ('_B07_' in x):
bands_dict['B07'] = R_dir+x
elif ('_B06_' in x):
bands_dict['B06'] = R_dir+x
elif ('_B05_' in x):
bands_dict['B05'] = R_dir+x
elif ('_B04_' in x):
bands_dict['B04'] = R_dir+x
elif ('_B03_' in x):
bands_dict['B03'] = R_dir+x
elif ('_B02_' in x):
bands_dict['B02'] = R_dir+x
else:
pass
return bands_dict
Selezione dell'area di studio (footprint)¶
- Identifico l'area di studio attraverso un
.geojsonnello specifico utilizziamo questo servizio online http://geojson.io/ - Cosa sono i file
.geojson? - Cos'è il formato wkt?
In [0]:
%%writefile ./geoj/venezia.geojson
{
"type": "FeatureCollection",
"features": [
{
"type": "Feature",
"properties": {},
"geometry": {
"type": "Polygon",
"coordinates": [
[
[
12.316017150878906,
45.46807675634549
],
[
12.284774780273438,
45.4324308667588
],
[
12.288551330566406,
45.340080499681875
],
[
12.344512939453125,
45.33284041773058
],
[
12.41455078125,
45.44471679159555
],
[
12.441329956054688,
45.47433654565603
],
[
12.316017150878906,
45.46807675634549
]
]
]
}
}
]
}
In [0]:
# Selezione del nostro .geojson file
GEOJSON_FILE = './geoj/venezia.geojson'
# Generazione del poligono in formato wkt
footprint = geojson_to_wkt(read_geojson(GEOJSON_FILE))
print(f"Stringa identificativa del poligono disegnato in formato wtk:\n{footprint}")
# Visualizzazione del poligono
bounds = folium.GeoJson(GEOJSON_FILE).get_bounds()
c_lat = (bounds[0][0] + bounds[1][0])/2
c_lon = (bounds[0][1] + bounds[1][1])/2
m = folium.Map([c_lat,c_lon], tiles='Stamen Terrain')
folium.GeoJson(GEOJSON_FILE).add_to(m)
m
Out[0]:
Sentinelsat API¶
- Iscriversi al sito ufficiale scihub https://scihub.copernicus.eu/dhus/#/home
- Premier sulle API
Credenziali¶
In [0]:
# Credenziali
USER = 'YOURUSR'
PSW = 'YOURPSW'
Inizializzazione API¶
In [0]:
# Inizializzazione oggetto API
api = SentinelAPI(USER, PSW, 'https://scihub.copernicus.eu/dhus')
# Parametri query
INIT_DATE = '20190101'
END_DATE = '20191115'
PLATFORM_NAME = 'Sentinel-2'
PROCESSING_LEVEL = 'Level-2A'
# Ricerca tramite area di studio
products = api.query(footprint,
date = (INIT_DATE, END_DATE),
platformname=PLATFORM_NAME,
processinglevel = PROCESSING_LEVEL,
cloudcoverpercentage = (0,90))
print(f"{len(products)} prodotti trovati.")
Conversione dei risultati in formato Geodataframe¶
- Cos'è un Geodataframe?
In [0]:
products_gdf = api.to_geodataframe(products)
products_gdf_sorted = products_gdf.sort_values(['cloudcoverpercentage'], ascending=[True])
print(products_gdf_sorted['cloudcoverpercentage'])
products_gdf_sorted
Out[0]:
Download prodotto ed estrazione .zip file¶
- Cosa c'è dentro alla folder?
In [0]:
# Selezione e Download del prodotto
product_uuid = products_gdf_sorted['uuid'][10]
product_uuid = '50ab5827-0e38-4a40-9a68-ade23b9e601c'
product_specs = products_gdf_sorted.loc[product_uuid]
print(f"{product_specs}\nPoligono del prodotto:")
product_specs.geometry
Out[0]:
Selezione dell'immagine a risoluzione desiderata¶
In [0]:
RESOLUTION = 10
if product_specs.filename in os.listdir("./data/"):
print("Prodotto già presente")
product_dir = F"./data/{product_specs.filename}"
else:
print(f"Download del prodotto {product_uuid} ...")
product_zip = api.download(product_uuid)
zip_name = product_zip['path']
print(f"Estrazione di {zip_name} ...")
with zipfile.ZipFile(zip_name, 'r') as zip_ref:
zip_ref.extractall("./data/")
os.remove(zip_name)
zip_name = f"./data/{zip_name[2:]}"
product_dir = zip_name.replace('.zip', '.SAFE')
granule_dir = f"{product_dir}/GRANULE/"
R_dir = f"{granule_dir}{os.listdir(granule_dir)[0]}/IMG_DATA/R{RESOLUTION}m/"
print(f"Cartella delle Risoluzioni {R_dir}")
Caricare le immagini (formato jpg2000)¶
- Che formato è?
- Che cosa sono le bande?
In [0]:
bands_dict = get_bands_files(R_dir)
print(f"Bande trovate")
bands_dict.keys()
Out[0]:
In [0]:
# Red
b4 = rio.open(bands_dict['B04'])
# Green
b3 = rio.open(bands_dict['B03'])
# Blue
b2 = rio.open(bands_dict['B02'])
# NIR
b8 = rio.open(bands_dict['B08'])
rim_info(b4)
In [0]:
fig, axs = plt.subplots(1,4, figsize=(30,30))
show(b2, ax=axs[0], title="B02", cmap="Blues")
show(b3, ax=axs[1], title="B03", cmap="Greens")
show(b4, ax=axs[2], title="B04", cmap="Reds")
show(b8, ax=axs[3], title="B08", cmap="plasma")
plt.show()
plt.close()
Creazione immagine .tiff RGB¶
In [0]:
with rio.open('./imgs/RGB.tiff','w',driver='Gtiff', width=b4.width, height=b4.height, count=3,crs=b4.crs,transform=b4.transform, dtype=b4.dtypes[0]) as rgb:
rgb.write(b2.read(1),3) # Red
rgb.write(b3.read(1),2) # Green
rgb.write(b4.read(1),1) # Blue
Clip dell'immagine all'area di studio¶
- Le immagini native sono troppo grandi, vanno clippate alla shape selezionata
In [0]:
# Dobbiamo proiettare le coordinate sul crs dell'immagine scaricata
footprint = gpd.read_file(GEOJSON_FILE)
footprint_proj = footprint.to_crs({'init': b4.crs})
# Clip dell'immagine
with rio.open("./imgs/RGB.tiff") as src:
out_image, out_transform = rio.mask.mask(src, footprint_proj.geometry ,crop=True)
out_meta = src.meta.copy()
out_meta.update({"driver": "GTiff",
"height": out_image.shape[1],
"width": out_image.shape[2],
"transform": out_transform})
with rio.open("./imgs/RGB_masked.tif", "w", **out_meta) as dest:
dest.write(out_image)
In [0]:
#@title Contrasto { run: "auto", vertical-output: true }
CONTRAST = 4 #@param {type:"slider", min:1, max:10, step:1}
with rio.open(r"./imgs/RGB_masked.tif") as msk:
# Correzione dati RGB
img = msk.read([1,2,3]).astype(np.float)
img[0] /= img[0].max()
img[1] /= img[1].max()
img[2] /= img[2].max()
fig, ax = plt.subplots(1, figsize=(12, 12))
show(np.clip((img*CONTRAST),0,1), title="Area di studio")
# Correzione canali
img_new = np.clip((img*CONTRAST*255),0,255).astype(np.uint8)
img_new = np.transpose(img_new, (1,2,0))
# Scrittura immagine
imageio.imsave("./imgs/risultato.jpg", img_new)
Sharpening Immagine¶
In [0]:
kernel_sharpening = np.array([[-1,-1,-1],
[-1, 9,-1],
[-1,-1,-1]])
sharpened = cv2.filter2D(img_new, -1, kernel_sharpening)
fig, ax = plt.subplots(1, figsize=(12, 12))
ax.imshow(sharpened)
plt.title("Immagine sharpened")
imageio.imsave("./imgs/risultato_sharpened.jpg", sharpened)
Calcolo Indici¶
In [0]:
nir = b8.read(1).astype(float)
green = b3.read(1).astype(float)
red = b4.read(1).astype(float)
ndvi = (nir-red)/(nir+red+1e-16)
In [0]:
plt.figure(figsize=(12,12))
show(ndvi+1, title="NDVI", cmap="Greens")
plt.close()
Creazione Video soglie¶
In [0]:
for i in range(10):
fig = plt.figure(frameon=False)
ax = plt.Axes(fig, [0., 0., 1., 1.])
ax.set_axis_off()
fig.add_axes(ax)
ax.imshow(ndvi>(i/10), cmap="Greens")
fig.savefig(f"./video/{i:03}.jpg", dpi=150)
plt.close()
!ffmpeg -y -loglevel panic -framerate 1 -i "./video/00%d.jpg" ./video/query_video.mp4
playvideo("./video/query_video.mp4")
Out[0]:
Clip NDVI al footprint¶
In [0]:
with rio.open('./imgs/ndvi.tiff','w',driver='Gtiff', width=b4.width, height=b4.height, count=1,crs=b4.crs,transform=b4.transform, dtype='float64') as ndvi_tiff:
ndvi_tiff.write(ndvi,1)
In [0]:
#@title NDVI Footprint { run: "auto", vertical-output: true }
SOGLIA = 0.9 #@param {type:"slider", min:0, max:1, step:0.1}
with rio.open("./imgs/ndvi.tiff") as src:
out_image, out_transform = rio.mask.mask(src, footprint_proj.geometry ,crop=True)
out_image = out_image[0,:,:]
out_image /= out_image.max()
out_image = np.clip(out_image,0,1)
plt.figure( figsize=(12,12))
plt.imshow(out_image>SOGLIA, cmap='Greens')
plt.title("NDVI footprint")
#plt.savefig("./imgs/ndvi_masked.jpg", dpi=300)
imageio.imwrite("./imgs/ndvi_masked.jpg", out_image)
In [0]:
# Plot heatmap NDVI
plt.figure(figsize=(15,15))
plt.imshow(out_image, cmap='Greens')
plt.title("Heatmap")
Out[0]:
Hands on LAB 👩💻¶
CVI¶
Ora tocca a voi...calcoliamo il Clorophyll Vegetation Index :)
- Implementare il calcolo dell'indice usando le bande
- Visualizzare l'heatmap dell'indice CVI
Calcolo Index:
$cvi = nir \cdot \frac{red}{green^{2}}$
Hints:
- Se non visualizzi il CVI put your Data Scientist hat :)
In [0]:
#-------------------------------Calcolo CVI-------------------------------------
cvi =
#-------------------------------------------------------------------------------
with rio.open('./imgs/cvi.tiff','w',driver='Gtiff', width=b4.width, height=b4.height, count=1,crs=b4.crs,transform=b4.transform, dtype='float64') as cvi_tiff:
cvi_tiff.write(cvi,1)
In [0]:
#@title CVI Footprint { run: "auto", vertical-output: true }
SOGLIA = 2.5 #@param {type:"slider", min:0, max:10, step:0.5}
with rio.open("./imgs/cvi.tiff") as src:
out_image, out_transform = rio.mask.mask(src, footprint_proj.geometry ,crop=True)
out_image = out_image[0,:,:]
plt.figure(figsize=(12,12))
plt.imshow(out_image>SOGLIA, cmap="Greens")
plt.title("CVI Footprint")
In [0]:
#-------------Riusciamo a visualizzare lheatmap delll'indice?-------------------
#-------------------------------------------------------------------------------
In [0]:
# Rilascio Risorse
# Attenzione se esegue questa cella, non è più possibile usare le bande
b2.close()
b3.close()
b4.close()
b8.close()
Now it's up to you :)
Thanks :)
Risorse¶
- Pangeo Project a Big Data Geoscience Community - https://pangeo.io/
- Database per ricerca di indici calcolabili dai satelliti https://www.indexdatabase.de/db/is.php
- Documentazione rasterio → manipolazioni immagini satellitari
- Documentazione geopandas → manipolazioni di dati geospaziali
- Documentazione folium → mappe interattive
- Documentazinoe sentinelsat → API per accedere ai dati dei sensori Sentinel