Propunem să prezentăm capabilitățiile catalogului de date creat recent în cadrul departamentului de Informatică al Universității de Vest din Timișoara. Catalogul expune date Sentinel-2 (în prezent) folosind cataloage STAC

Stocăm date Sentinel-2 la nivelul de procesare L1C din 2017 până în prezent. Date care, acoperă întreaga suprafață a României. În cadrul acestui articol, vom ilustra folosirea bibliotecii stackstac pentru transformarea colecției STAC într-un dataframe xarray care va fi procesat folosind dask pentru generarea unor GIF-uri ce arată starea unor orașe din România de alungul aniilor. Lucrarea prezentata aici se bazeaza pe un exemplu stackstac

Abordarea pentru generarea acestor GIF-uri este destul de simplă:

1. Generăm un areal cu ajutorul codurilor SIRUTA folosindu-ne de către un serviciu pus la dispoziție de geo-spatial.org ce conține limitele administrativ teritoriale ale României.
2. Interogăm catalogul de date folosind stackstac și selectăm benziile RGB pentru Sentinel-2 (B02, B03 respectiv B04) și convertim aceste date într-un array xarray.
3. Filtrăm scenele pe baza acoperirii acestora cu nori (metadatele acestea fiind disponibile în catalogul STAC). Pentru exemplele de mai jos, scenele au fost filtrate astfel încât oricare dintre acestea sa aiba în jur de 1% acoperire cu nori.
4. Redimensionarea array-ului xarray de-alungul dimensiunii de timp la fiecare 10 zile.
5. Decuparea array-ului cu ajutorului unui buffer (exprimat în metrii) față de punctul din centrul arealului.
6. Aplicarea unei corecții gamma pentru crearea unui compozit True Color Image (TCI), folosind o valoare de 2.2 pentru parametrul gamma (Rezultatele prezentate într-o lucrare publicata de compania Sinergise arată faptul că aceasta este o valoare ce funcționează bine pentru date Sentinel-2)
7. Folosirea biblioteciilor dask și geogif pentru randarea GIF-urilor.

Cum abordam problema în Python?

începem prin interogarea catalogului STAC:

1
2
3
4
5
6
7
8
9

url = 'https://registry.sage.uvt.ro/geoserver/ogc/stac?f=application%2Fgeo%2Bjson'

catalog = pystac_client.Client.open(url)

result = catalog.search(
    bbox=bounds_wgs84,
    collections=["S2L1C"],
    datetime='2017-01-12T00:00:00Z/2022-05-31T12:31:12Z',
)

Folosim stackstac.stack() pentru crearea array-ului:

1

stack = stackstac.stack(items, epsg=32635, resolution=10, assets=["B04", "B03", "B02"], errors_as_nodata=Exception(".*"))

Acesta este cubul de date creat:

Putem apoi filtra scene pe baza acoperirii acestora cu nori

1
2
3

stack = stack[stack["eo:cloud_cover"] < 1]

composites = stack.resample(time="10D").median("time")

Convertim apoi coordonatele lat-lon ale punctului din centrul zonei administrative teritoriale în proiecția UTM cu ajutorul bibliotecii proj și decupăm arealul de interes:

1
2
3
4

x_utm, y_utm = pyproj.Proj(stack.crs)(lat, lon)
buffer = 5000 # meters

aoi = composites.loc[..., y_utm+buffer:y_utm-buffer, x_utm-buffer:x_utm+buffer]

După care putem folosi dask.array.power() pentru a aplica corecția gamma pe cubul de date:

1
2

def gamma_correction(xarr, gamma):
    return da.power(xarr, 1/gamma)

Ca un ultim pas, ne folosim de dask și geogif pentru producerea GIF-urilor

1
2
3

with dask.diagnostics.ProgressBar():
    gif_img = geogif.dgif(aoi, fps=7, date_color=(255,255,255), date_bg=(0,10,0)).compute()
    dask.utils.format_bytes(len(gif_img.data))

Acestea sunt rezultatele:

Doresc sa îi mulțumesc colegului meu, Marian Neagul pentru sprijinul si îndrumarea oferită, organizației geo-spatial.org pentru punerea la dispoziție a limitelor administrativ-teritoriale (și a unor multitudini de alte date) dar și Agenției Spațiale Române pentru punerea la dispoziție a unui mediu Jupyter Hub.