Věděli jste, že nám patří celá letka družic, které nepřetržitě snímkují celý svět? A protože jsou opravdu naše, jejich čerstvé záběry si teď hned můžete zdarma stáhnout a libovolně využít.
Řeč je samozřejmě o programu dálkového průzkumu Země Evropské unie Copernicus a jeho družicových misích Sentinel.
Nad hlavami nám poletují už od roku 2014, a pokud nevyschnou evropské penězovody, v následujících letech přibydou další. V závěru desetiletí snad i mise Sentinel-12 vyzbrojená radarem SAR, který vyfotí zemský povrch dokonce i při silné oblačnosti, jeho paprsky totiž projdou skrze zamračenou oblohu.
Windy a vrstva s koncentrací NO₂ (především spalovací motory), kterou aplikace čerpá právě z obrovských databází programu Copernicus
Za zmínku stojí také mise Sentinel-5, kterou velmi dobře znají fanoušci české meteorologické mapy Windy. Družice je vyzbrojená hromadou čidel pro dálkové měření kvality ovzduší a Windy její data zobrazuje třeba na vrstvě NO₂ (emise spalovacích motorů).
Sentinel-2 snímá planetu v rozlišení 10 metrů/pixel a přinejmenším 1× týdně
My se dnes ale spojíme ještě s jednou misí – nejspíše tou zdaleka nejatraktivnější –, budeme totiž stahovat snímky z družic rodiny Sentinel-2 vyzbrojených zařízením MSI (MultiSpectral Instrument).
Evropská družice z rodiny Sentinel-2
MSI je speciální fotoaparát, který snímá zemský povrch v širokém spektru vlnových délek včetně přirozeného světla a družice vyfotí prakticky libovolné místo v Evropě zhruba každý druhý-třetí den při rozlišení až 10 metrů/pixel.
Snímací MultiSpectral Instrument na palubě družice Sentinel-2
Stručně řečeno, díky Sentinel-2 si můžete prohlédnout čerstvý družicový snímek vašeho města, který bude zpravidla nejvýše několik málo dnů starý. A přesně to si dneska vyzkoušíme a zautomatizujeme v Pythonu!
80+ PB dat z družic Sentinel
Copernicus za ty roky nasbíral více než 80 petabajtů dat, ke kterým můžeme zdarma přistupovat a zpracovávat je pro libovolné účely. A protože je toho opravdu hodně, během doby se zrodilo hned několik prohlížečů a zrcadel. Oficiálním okýnkem do světa družic Sentinel je ale Copernicus Browser.
Copernics Browser a pohled na Prahu z družice Sentinel-2 při přeletu v pátek 6.9. 2024. Stejný snímek si dnes stáhneme i pomocí Pythonu a oficiálního API pro vývojáře
I ten je samozřejmě bezplatný, pro plnohodnotné použití se ale musíte zaregistrovat. Poté se podobně jako v jakémkoliv jiném mapovém portálu posunete třeba nad Prahu a v postranním panelu dohledáte nejčerstvější snímek.
Služba umí snímky také filtrovat, zobrazovat různé zachycené vlnové délky, ve kterých lépe vynikne vodstvo, zdravá zemědělská půda nebo třeba sněhová pokrývka, a také zvládne skládat časosběrná videa a provádět nejrůznější měření.
Jak se dostat k datům z družic úplně zadarmo
K surovým datům nakonec můžeme přistupovat i skrze nespočet aplikačních rozhraní pro vývojáře. I zde platí, že jich je celá řada a stejně tak existuje hned několik zrcadel archivu, z nichž některá jsou komerční. To je případ zejména (skvělé) slovinské služby Sentinel Hub, která má k dispozici i knihovnu pro snadnou práci v Pythonu.
My to všechno ale chceme opravdu zadarmo, a tak se dnes spojíme přímo s originálem, který se jmenuje CDSE – Copernicus Data Space Ecosystem. Je to obrovský archiv prakticky všech dat z mnohaletých misí a i k němu můžeme přistupovat skrze mnohá rozhraní.
Catalog API pro hledání a Process API pro tvorbu snímků
Jedno z rozhraní se jmenuje (také) Sentinel Hub API, je totiž kompatibilní s výše zmíněnou cizí službou, ale protože zůstaneme na serverech Koperníka, bude to zadarmo. Ano, je to matoucí, ale co v Evropě není, viďte…
Sentinel Hub API používá běžné HTTP GET/POST a několik dotazů. My si dnes vyzkoušíme:
I když jsou obě API zdarma, opravdu si musíte vytvořit účet, budete totiž potřebovat přístupový klíč, aby vás mohl Copernicus případně dočasně zablokovat, pokud přesáhnete kvótu při masivním stahování celé planety.
Copernicus Request Builder je skvělý interaktivní pomocník při studiu a zkoušení API
Práci se všemi API z rodiny Sentinel Hub si můžete vyzkoušet na interaktivním webu Copernicus Request Builder, kde si celý HTTP požadavek postavíte na míru a generátor poté ukáže kód pro několik běhových prostředí. Třeba pro Python s využitím knihovny Request a pro textový HTTP klient cURL.
Autentizace API pomocí OAuth
Tak, dost formalit, pojďme si to vyzkoušet v Pythonu. Výše zmíněný Request Builder po přihlášení sám doplnil API klíč/token, ten je ale dočasný, takže se při každém požadavku autentizujeme pomocí techniky OAuth. Nebude to nic složitého, pomůže nám totiž knihovna Requests-OAuthlib pro Python, která je nadstavbou pro knihovnu HTTP komunikace Requests.
Obě můžeme nainstalovat třeba pomocí správce balíčků pro Python pip:
pip install requests requests-oauthlib
V dalším kroku si musíme vytvořit údaje pro OAuth, což provedeme v nastavení účtu Copernicus. Vpravo v sekci OAuth clients vytvořte nového jako na obrázku níže, a pokud ho nebudete chtít obnovovat stále dokola, zaškrtněte volbu „Never expire.“
V nastavení účtu na púortálu Copernicus si vytvoříme OAuth klienta s přihlašovacími údaji
Formulář vám poté vygeneruje dva ověřovací údaje Client ID a Client secret, které použijeme pro autentizaci a vytvoření API klíče/tokenu.
Ahoj API, máš nějaké přelety nad Prahou z posledních dnů?
Fajn, máme přihlašovací údaje pro komunikaci, takže si pojďme vyzkoušet katalogové API. Dejme tomu, že chceme zjistit, jestli jsou na serveru nějaké snímky z mise Sentinel-2 L2A (záběry v přirozených barvách) s těmito parametry:
- Oblast: Praha
- Časový rozsah: 1. září 2024 až 7. září 2024 (včetně)
Copernicus samozřejmě netuší, kde je Praha, ve vyhledávání proto používáme souřadnice obdélníku (bbox), nebo složitějších polygonů, které zahrnují oblast našeho zájmu.
Pokud neznáte souřadnice, s výběrem obdélníku (bbox) pomůže web bboxfinder.com
Pomůže nám třeba web bboxfinder.com, kde si od ruky označíme okolí Prahy a v zápatí se nám zobrazí souřadnice obdélníku (jeho protilehlé rohy):
- 13.960876,49.841525,14.988098,50.296358
Časový rozsah se pak do API posílá ve standardním textovém formátu a v časové zóně na nultém poledníku:
- 2024-09-01T00:00:00Z/2024-09-07T23:59:59Z
Katalogové API prozradí, kdy Sentinel-2 přeletěl nad Prahou a jestli bylo hezké počasí
Všechny tyto údaje vložíme do struktury JSON a tu pošleme skrze katalogové API na server. Odpovědí bude opět sáhodlouhý JSON s nalezenými snímky, ve kterých nechybí ani velmi praktický procentuální údaj o míře oblačnosti. Díky tomu víme, jestli bude v nalezených datech vůbec vidět povrch, anebo jen zakaboněná zářijová obloha.
Katalog dat z družice Sentinel-2 L2A našel dva přelety nad Prahou ze středy a pátku. Oba proběhly při velmi nízké oblačnosti, takže to budou jistě velmi hezké záběry!
Kompletní kód pro vyhledání dat v Pythonu pomocí Catalog API a s autentizací pomocí OAuth, aniž bychom museli ručně vyplňovat dočasný token, bude vypadat takto:
from oauthlib.oauth2 import BackendApplicationClient
from requests_oauthlib import OAuth2Session
# Prihlasovaci udaje pro OAuth
# Vytvor zde: shapps.dataspace.copernicus.eu/dashboard/#/account/settings
client_id = ""
client_secret = ""
if len(client_id) == 0 or len(client_secret) == 0:
print("Vypln client_id a client_secret pro OAuth autentizaci!")
print("Zaloz si bezplatny ucet na Copernicus a bez sem: https://shapps.dataspace.copernicus.eu/dashboard/#/account/settings")
exit()
# Parametry API dotazu
parametry_dotazu = {
"bbox": [13.960876,49.841525,14.988098,50.296358],
"datetime": "2024-09-01T00:00:00Z/2024-09-07T23:59:59Z",
"collections": ["sentinel-2-l2a"],
"limit": 5
}
# OAuth autentizace
client = BackendApplicationClient(client_id=client_id)
oauth = OAuth2Session(client=client)
token = oauth.fetch_token(token_url="https://identity.dataspace.copernicus.eu/auth/realms/CDSE/protocol/openid-connect/token", client_secret=client_secret, include_client_id=True)
# Dotaz na Catalog API
odpoved = oauth.post("https://sh.dataspace.copernicus.eu/api/v1/catalog/1.0.0/search", json = parametry_dotazu)
if odpoved.status_code == 200:
data = odpoved.json()
if "features" in data:
print(f"Nalezeno zaznamu: {len(data['features'])}:")
for i,feature in enumerate(data["features"]):
print(f"{i}. zaznam:")
print(f"\tDatum: {feature['properties']['datetime']}")
print(f"\tPlatforma: {feature['properties']['platform']} (Zarizeni: {','.join(feature['properties']['instruments'])})")
print(f"\tZakryti mraky: {feature['properties']['eo:cloud_cover']} %")
print("")
else:
print(f"Chybicka se vloudila!\nHTTP kod: {odpoved.status_code}\n{odpoved.text}")
A teď to pojďme stáhnout jako JPEG obrázek!
Výborně, takže už umíme vyhledávat v katalogu a teď si ukážeme, jak se vlastně snímky stahují. Použijeme k tomu druhý typ dotazu jménem Process API. Jak už název napovídá, snímky vlastně nebudeme stahovat z nějaké mezipaměti, jak to funguje u běžných mapových portálů včetně Mapy.cz a Mapy Google, ale budeme je přímo generovat ze surových dat.
Čerstvý záběr vodního díla Nové Mlýny, který jsme si vygenerovali a stáhnuli pomocí našeho programu v Pythonu, jehož kód najdete v závěru článku
Na serverech Copernicus totiž nejsou uložené jakési hotové JPEG obrázky ve formě nařezaných dlaždic, ale opravdu jen zdrojové informace v různých kanálech zaznamenaného spektra.
Samotný dotaz tedy bude v principu hodně podobný jako u katalogového API, nicméně jeho součástí bude ještě Evalscript! Je to vlastně vcelku jednoduchý Javascript s podprogramem, který sděluje generátoru na serveru, jak má zpracovat surová data.
Jak už jsme si řekli výše a potvrdil to i výsledek katalogového API, družice Sentinel-2 jsou vyzbrojené zařízením MSI (MultiSpectral Instrument), které snímá povrch multispektrálně v různých vlnových délkách – kanálech – a ty jsou v Evalscrpitu řádně očíslované.
Vlnové délky, se kterými pracuje Sentinel-2 a jejich identifikátory pro evalskript
Každý slouží k něčemu jinému. Jedny zobrazí oblast v přirozených barvách, jiné umocní třeba vlhkost v půdě, což je klíčové pro zemědělce atp. Seznam všech kanálů se stručným vysvětlením, k čemu se hodí, najdete třeba zde, anebo více lidštěji také v této nápovědě.
Skládáme červenou, zelenou a modrou
Protože nás zajímají snímky v přirozených barvách RGB a nebudeme se pouštět do složitějších environmentálních analýz, v evalskriptu jen složíme tyto kanály:
- B02 (492 nm, modrá barva)
- B03 (560 nm, zelená barva)
- B04 (665 nm, červená barva)
Tyto kanály mají zároveň nejvyšší rozlišení 10 metrů/pixel, zatímco třeba u takového kanálu B01 (20 nm, citlivý na přítomnost aerosolu) je to pouze 60 metrů/pixel.
Surová data jsou zpravidla velmi tmavá – podexponovaná, takže je můžeme v evalskriptu ještě zesílit. Ale pozor, pokud budou v záběru třeba plechové střechy, sníh a silná oblačnost, znásobením můžeme vytvořit silný přepal.
Zároveň zpracováváním dat spotřebováváme kvótu, takže dodatečné vylepšení snímku se asi hodí provádět až po stažení. V každém případě, evalskript pro získání snímku v přirozených kanálech RGB a s dvojnásobným zesílením expozice bude vypadat takto:
//VERSION=3
function setup() {
return {
input: ["B02", "B03", "B04"],
output: { bands: 3 },
}
}
function evaluatePixel(sample) {
let expozice = 2;
return [sample.B04 * expozice, sample.B03 * expozice, sample.B02 * expozice];
}
Evalskript se posílá jako další parametr dotazu, ve kterém tentokrát oproti katalogovému dotazu vyplníme i formát výstupních dat. Může to být JPEG, ale i PNG a GeoTIFF pro další zpracování v nástrojích pro GIS a kartografii. A také JSON. Do toho se už ale pouštět nebudeme.
Na výstupu zároveň můžeme stanovit, jak má být obrázek velký, přičemž API dovolí vytvořit nejvýše bitmapu o šířce 2 500 pixelů.
Kompletní program pro stahování snímků
Procesní API nabízí hromadu dalších parametrů, které se nám už do článku nevejdou, všechny do jednoho si ale můžete opět vyzkoušet v testovacím Copernicus Request Builder.
Na závěr si proto už rovnou ukážeme hotovou aplikaci pro stahování nejmladšího snímku pro zvolenou oblast. Jako kritérium výběru použijeme časový filtr pro posledních sedm dnů, abychom měli rezervu pro všechny případy.
Družice Sentinel-2 sice formálně snímají s periodou až 5 dnů, jak už jsme si ale řekli výše a ukázali i v katalogovém API, v praxi je v Evropě frekvence spíše vyšší a snímek Prahy nebo Brna se na serveru objeví zpravidla každé 2-3 dny.
Program budeme ovládat pomocí několika přepínačů/parametrů:
- -S, --soubor (název výstupního souboru, třeba praha.jpg)
- -V, --vyska (výška obrázku v px, šířku dopočítáme automaticky)
- -B, --box (souřadnice výběrového obdélníku)
- -F, --format (formát výstupních dat)
- -K, --kvalita (kvalita obrázku, 1-100)
- -E, --evalscript (cesta k vlastnímu volitelnému evalskriptu)
- -J, --jas (násobek zvýšení jasu na serveru, třeba 4)
- -T, --ukaztoken (pouze se přihlásíme a vypíšeme token)
Pojďme si to vyzkoušet naostro. Dejme tomu, že chceme získat JPEG obrázek Prahy jako soubor praha.jpg podle výběrového obdélníku, který jsme použili už o katalogového API. Zvýšíme jas obrázku 5×, protože surová data jsou při nízké oblačnosti silně podexponovaná, a výšku nastavíme na 1000 pixelů.
Šířka se dopočítá automaticky podle poměru stran výběrového obdélníku a výška případně sníží, pokud by byla šířka příliš velká.
Spouštíme náš program a za pár sekund už máme čerstvý snímek Prahy a okolí
Pokud pojmenujeme program jako sentinel2.py, pak bude vypadat povel ke stažení obrázku takto:
python .\sentinel2.py --soubor praha.jpg --box 13.960876,49.841525,14.988098,50.296358 --vyska 1000 --format image/jpeg --jas 4
Opět připomínám, že pro správný běh programu musíte mít nainstalované běhové prostředí Pythonu a nadstandardní knihovny Requests a Requests-OAuthlib!
Nebudou-li servery zrovna přetížené, dohledávání nejčerstvějších dat a generování obrázku zabere nejvýše jednotky sekund a za okamžik už vám na počítači přistane kýžený obrázek.
Praha a okolí na snímku z 6. září 2024. Záběr je tedy starý pouhých 24 hodin!
Pro sorvnání se podívejte, jak by vypadal surový snímek bez navýšení jasu (--jas 1) a naopak s ještě více zesílenou expozicí (--jas 6).
Hodnotu vždy vyzkoušejte nad konkrétním území, odrazy slunce z plechových střech totiž mohou způsobit ohromné přepaly.
Surové kanály B02, B03 a B04 v původních hodnotách a s šestinásobným zesílením
Skript můžete volat třeba jednou týdně
Skript poté můžete automaticky spouštět třeba 1-2 do týdne a sledovat, jak se mění konkrétní území bez ruční kontroly v Copernicus Browser.
A protože si můžete libovolně upravit i procesní evalskript – koukněte na tyto příklady –, snadno si vyrobíte třeba detektor, který stáhne obrázky v umělých barvách třeba pro měření kvality vegetace, sucha na zemědělské půdě atp.
Kompletní zdrojový kód
Tak jo, máme hotovo a lovení družicových snímků zdar! Na úplný závěr samozřejmě nesmí chybět to nejdůležitější, tedy kompletní kód naší stahovací miniaplikace v Pythonu.
Skript najdete i na GitHubu našeho seriálu o programování elektroniky.
# Pro OAuth se pouziva https://oauthlib.readthedocs.io/en/latest/installation.html
# Nainstaluj smudlo!
from oauthlib.oauth2 import BackendApplicationClient
from requests_oauthlib import OAuth2Session
from datetime import datetime, timedelta, timezone
import math, argparse
# Nastaveni API OAuth (musis se zaregistrovat)
# https://shapps.dataspace.copernicus.eu/dashboard/#/account/settings
client_id = ""
client_secret = ""
if len(client_id) == 0 or len(client_secret) == 0:
print("Vypln client_id a client_secret pro OAuth autentizaci!")
print("Zaloz si bezplatny ucet na Copernicus a bez sem: https://shapps.dataspace.copernicus.eu/dashboard/#/account/settings")
exit()
# Vychozi hodnoty
VYSKA = 512 # Vychozi vyska obrazku
SIRKA_MAX = 2500 # Maximalni povolena sirka obrazku na strane API
FORMAT = "image/jpeg" # Format vystupnich dat
SOUBOR_NAZEV = "snimek.jpg" # Vystupni soubor
KVALITA = 90 # JPEG kvalita
BOX = "13.960876,49.841525,14.988098,50.296358" # GPS ouradnice ctverce Prahy
# Evalscript je specialni procesni Javascript, ktery se spousti nad daty na serverech Copernicus
# Chceme druzicova data v pravych barvach (kanaly B02, B03, B04)
# Zesilime expozici, surova data jsou totiz hodne tmava
EVALSCRIPT = """
//VERSION=3
function setup() {
return {
input: ["B02", "B03", "B04"],
output: { bands: 3 },
}
}
function evaluatePixel(sample) {
let expozice = 2;
return [sample.B04 * expozice, sample.B03 * expozice, sample.B02 * expozice];
}
"""
parser = argparse.ArgumentParser()
parser.add_argument("-S", "--soubor", default=SOUBOR_NAZEV, type=str, help = "Nazev souboru (praha.jpg apod.)")
parser.add_argument("-V", "--vyska", default=VYSKA, type=int, help = "Vyska obrazku v pixelech (500, 1000 apod.)")
parser.add_argument("-B", "--box", default=BOX, type=str, help = "Vyberovy obedlnik v souradnicich WGS84 (lng0,lat0,lng1,lat1)")
parser.add_argument("-F", "--format", default=FORMAT, type=str, help = "Format obrazku (image/jpeg, image/png)")
parser.add_argument("-K", "--kvalita", default=KVALITA, type=int, help = "Kvalita obrazku (0-100)")
parser.add_argument("-E", "--evalscript", default="", type=str, help = "Soubor vlastniho evalskriptu (evalscript.js)")
parser.add_argument("-J", "--jas", default=2, type=float, help = "Jas obrazku, pokud se pouzije vestaveny evalskript (0-XXX)")
parser.add_argument("-T", "--ukaztoken", action="store_true", help = "Pouze vygeneruje docasny token/API klic (lze pouzit pro Bearer autentizaci bez OAuth)")
argumenty = parser.parse_args()
# Uprava parametru podle argumentu skriptu
VYSKA = argumenty.vyska
FORMAT = argumenty.format
SOUBOR_NAZEV = argumenty.soubor
KVALITA = argumenty.kvalita
BOX = [float(souradnice) for souradnice in argumenty.box.split(",")]
EVALSCRIPT = EVALSCRIPT.replace("let expozice = 2;", f"let expozice = {argumenty.jas};")
# Pokud chceme pouzit vlastni evalskript ze souboru, prepiseme ten vychozi
if len(argumenty.evalscript) > 1:
try:
with open(argumenty.evalscript, "r") as f:
EVALSCRIPT = f.read()
except:
print(f"Nemohu přečíst evalscript {argumenty.evalscript}. Použiji výchozí")
# Vytvorime OAuth sezeni
client = BackendApplicationClient(client_id=client_id)
oauth = OAuth2Session(client=client)
token = oauth.fetch_token(token_url='https://identity.dataspace.copernicus.eu/auth/realms/CDSE/protocol/openid-connect/token', client_secret=client_secret, include_client_id=True)
# Pokud chceme jen vypsat token/API key pro Bearer autentizaci
# Vypiseme autentizacni hodnoty a ukoncime program
if argumenty.ukaztoken:
for key in token:
print(f"{key}: {token[key]}")
exit()
# Vypocet pomeru stran pro obrazek v pixelech ze zemepisnych souradnic
# Potencialne nepresne, ale jde nam o obrazek pro bezne pouziti, ne geodezii a GIS
rozdil_sirky = BOX[3] - BOX[1]
rozdil_delky = BOX[2] - BOX[0]
prumerna_sirka = (BOX[1] + BOX[3]) / 2
vzdalenost_sirky = rozdil_sirky * 111
vzdalenost_delky = rozdil_delky * 111 * math.cos(math.radians(prumerna_sirka))
pomer_stran = vzdalenost_delky / vzdalenost_sirky
# Vypocet sirky vystupniho obrazku podle nastavene vysky a poměru stran
sirka = int(VYSKA * pomer_stran)
vyska = VYSKA
# API umi vytvorit obrazek o sirce nejvyse 2500 pixelu
# Takze pokud jsme vypocitali vetsi sirku, zmensime vysku
if sirka > SIRKA_MAX:
sirka = SIRKA_MAX
vyska = int(sirka / pomer_stran)
print(f"Rozmery hotoveho obrazku: {sirka}x{vyska} pixelu")
# Druzice nad mistem/boxem v Evrope sice leti zpravidla kazde 2-3 dny,
# ale prodleva muze byt i delsi, a tak stanovime vyhledavaci rozsah na 7 dnu
datum_konec = datetime.now(timezone.utc)
datum_zacatek = (datum_konec - timedelta(days=7)).isoformat()
datum_konec = datum_konec.isoformat()
# Priprava API dotazu ve forme JSON
# Budeme chit snimek z druzice Sentinel-2 L2A
# Pro dany box, pixelovou velikost, format atp.
# Pro dalsi mozne parametry kouknete na https://shapps.dataspace.copernicus.eu/requests-builder/
parametry_dotazu = {
"input": {
"bounds": {
"bbox": BOX
},
"data": [
{
"type": "sentinel-2-l2a",
"dataFilter": {
"timeRange": {
"from": datum_zacatek,
"to": datum_konec,
}
},
"processing": {
"upsampling": "BILINEAR",
"downsampling": "BILINEAR"
},
}
]
},
"output": {
"width": sirka,
"height": vyska,
"responses": [
{
"identifier": "default",
"format": {
"type": FORMAT,
"quality": KVALITA,
}
}
]
},
"evalscript": EVALSCRIPT
}
# Provedeme autentizovany HTTP dotaz skrze OAuth a server by mel odpovedet rovnou bajty obrazku
# V tomto API dotazu se (asi) nedozvime, z jakeho presneho casu ve vyberu obrazek pochazi
odpoved = oauth.post("https://sh.dataspace.copernicus.eu/api/v1/process", json = parametry_dotazu)
if odpoved.status_code == 200:
print(f"Ukladam obrazek {SOUBOR_NAZEV}")
with open(SOUBOR_NAZEV, "wb") as file:
file.write(odpoved.content)
else:
print(f"Chybicka se vloudila!\nHTTP kod: {odpoved.status_code}\n{odpoved.text}")