ETa Mapper
Earth Engine
Statisztika Statistics
GeoTIFF

Napi ETa-térképek Sentinel‑2 és ERA5‑Land alapján Daily ETa maps from Sentinel-2 and ERA5-Land

  1. Jelöld ki a területet (alapból a Kis‑Balaton) Choose the area (Kis‑Balaton by default)
  2. Állítsd be az időszakot Set the period
  3. Indítsd a számítást Start the calculation
Poligon-statisztika Polygon statistics

ETa a megrajzolt poligonok alatt, minden napra ETa under the drawn polygons, every date

Számolás után rajzolj poligont vagy téglalapot a térkép fölötti Statisztika gombokkal. Itt minden napra megjelenik alatta az ETa átlaga és mediánja. After a calculation, draw a polygon or rectangle with the Statistics buttons above the map. The mean and median ETa under it appear here for every date.

Módszertan Method

Műholdképből párolgástérkép – lépésről lépésre From satellite image to evaporation map – step by step

Ezt az appot arra írtam, hogy a GEE_code.txt számítását bárki lefuttathassa egy kattintással, tetszőleges területre és időszakra. Közben átnéztem a scriptet, és hat ponton javítottam rajta. Itt végigvezetlek azon, mit csinál az app, és hol és miért nyúltam az eredeti kódhoz. Az ábrák interaktívak – kattints, húzz, kapcsolgass. I wrote this app so that anyone can run the calculation of GEE_code.txt with one click, for any area and period. Along the way I reviewed the script and corrected it in six places. Here I walk you through what the app does, and where and why I changed the original code. The figures are interactive – click, drag, toggle.

8 lépés steps
6 javítás corrections
~2,5× ~2.5× túlbecslés az eredetiben overestimate in the original
100% pixelegyezés az eredeti TIF-ekkel pixel match with the original TIFs
01 lépés step

Kijelölöd a területet és az időszakot You choose the area and the period

Minden egy téglalappal vagy poligonnal kezdődik. Alapból a GEE_code.txt Kis-Balaton téglalapját használom (17,00–17,40° K, 46,50–46,80° É), de rajzolhatsz sajátot, vagy feltölthetsz egy GeoJSON-t, KML-t, GeoPackage-et vagy shapefile-t.

Az időszakon belül megkeresem az összes napot, amikor a Sentinel-2 műholdak (A, B, C) elhaladtak a terület felett – Magyarországon ez nagyjából 2–3 naponta van. A felhőborítottság-határ kiszűri azokat a felvételeket, amelyek túlságosan felhősek. Húzd a csúszkát az ábrán, és nézd, hány nap marad!

Everything starts with a rectangle or a polygon. By default I use the Kis-Balaton rectangle of GEE_code.txt (17.00–17.40° E, 46.50–46.80° N), but you can draw your own or upload a GeoJSON, KML, GeoPackage or shapefile.

Within the period I look for every day on which a Sentinel-2 satellite (A, B, C) passed over the area – in Hungary roughly every 2–3 days. The cloud limit drops scenes that are too cloudy. Drag the slider in the figure and see how many days remain!

02 lépés step

Napi mozaikot készítek, és kimaszkolom a felhőket I build a daily mosaic and mask the clouds

Egy Sentinel-2 felvétel 110 × 110 km-es csempe. A Kis-Balaton éppen két csempe határán fekszik, így egy napon gyakran két darabban érkezik a kép. Az eredeti script egyenként, az első 10 felvételt dolgozta fel; én naponta egy mozaikba rakom őket, így egy dátum = egy, a teljes területet lefedő kép.

A javított módban a SCL (scene classification) sáv alapján pixelenként kidobom a felhőt, a felhőárnyékot és a cirruszt. Az eredeti script csak a felvétel egészére vonatkozó felhő-százalékot nézte, így egy 9%-ban felhős kép felhős pixelei is ETa-értéket kaptak.

A Sentinel-2 scene is a 110 × 110 km tile. Kis-Balaton lies right on the border of two tiles, so on one day the image often comes in two pieces. The original script processed the first 10 scenes one by one; I put them into one mosaic per day, so one date = one image covering the whole area.

In the corrected mode I remove cloud, cloud shadow and cirrus pixel by pixel with the SCL (scene classification) band. The original script only looked at the cloud percentage of the whole scene, so the cloudy pixels of a 9% cloudy scene still got an ETa value.

03 lépés step

Kiszámolom a növényzet zöldségét: NDVI I measure how green the vegetation is: NDVI

Az NDVI a közeli infravörös (B8) és a vörös (B4) visszaverődés normalizált különbsége: NDVI = (B8 − B4) / (B8 + B4). Az egészséges növény a vöröset elnyeli, az infravöröst visszaveri, ezért értéke 0,6–0,9; a csupasz talaj 0,1–0,2 körül van, a víz 0 alatt.

Egy csapda: 2022. január 25. óta az ESA a nyers értékekhez +1000-et ad. A régi COPERNICUS/S2_SR gyűjtemény ezt nem vonja le, így az NDVI lefelé torzul. Kapcsold be az ábrán az eltolást, és nézd, mennyit csökken! Ezért a javított módban a S2_SR_HARMONIZED gyűjteményt használom.

NDVI is the normalised difference of near-infrared (B8) and red (B4) reflectance: NDVI = (B8 − B4) / (B8 + B4). Healthy plants absorb red and reflect infrared, giving 0.6–0.9; bare soil is around 0.1–0.2 and water is below 0.

A trap: since 25 January 2022 ESA adds +1000 to the raw values. The old COPERNICUS/S2_SR collection does not remove it, so NDVI is biased low. Switch the offset on in the figure and see how much it drops! That is why the corrected mode uses S2_SR_HARMONIZED.

04 lépés step

Az NDVI-ből növényi együtthatót (Kc) rendelek I turn NDVI into a crop coefficient (Kc)

A Kc megmondja, egy felszín a referencia-gyephez képest mennyit párologtat. A scriptben ez egy lépcsős táblázat az NDVI osztályai szerint – ezt változatlanul hagytam, és a config.R-ben átírható.

Húzd az NDVI-csúszkát: a pont végigfut a lépcsőn. Figyeld a 0,1–0,2-es osztályt! A scriptben ez „víz” megjegyzéssel 1,05-ös Kc-t kap, pedig a víz NDVI-je jellemzően 0 alatti, a 0,1–0,2 inkább ritka növényzet vagy nedves talaj. Ezt nem javítottam, mert szakmai döntés, de érdemes átgondolni.

Kc tells how much a surface evaporates compared with a reference grass. In the script it is a step table of NDVI classes – I kept it unchanged, and it can be edited in config.R.

Drag the NDVI slider: the dot runs along the steps. Watch the 0.1–0.2 class! The script gives it Kc = 1.05 with the comment “water”, although water usually has NDVI below 0; 0.1–0.2 is more like sparse vegetation or wet soil. I did not change this because it is a scientific choice, but it is worth reconsidering.

05 lépés step

Kiszámolom a referencia-párolgást (ET0) – itt volt a legtöbb hiba I compute reference evapotranspiration (ET0) – where most of the errors were

Az ET0 egy jól vízellátott, rövid füves felszín párolgása, a FAO-56 Penman–Monteith egyenlettel (Allen et al., 1998). Az ERA5-Land napi hőmérsékletéből, harmatpontjából, sugárzásából, szeléből és légnyomásából számolom.

Az ábra egy élő kalkulátor: állítsd be egy nyári nap időjárását, és kapcsold be egyenként a javításaimat. Látod, ahogy az eredeti script ~14 mm/napos értéke a reális ~5–6 mm/napra csökken. A legnagyobb tétel a hiányzó 0,408-as szorzó: ez váltja át a sugárzást MJ m−2-ből párolgó vízmilliméterré, nélküle a sugárzási tag 2,45-szörös.

ET0 is the evapotranspiration of a well-watered short grass, from the FAO-56 Penman–Monteith equation (Allen et al., 1998). I compute it from the daily temperature, dew point, radiation, wind and pressure of ERA5-Land.

The figure is a live calculator: set the weather of a summer day and switch my corrections on one by one. You can watch the ~14 mm/day of the original script drop to a realistic ~5–6 mm/day. The biggest item is the missing 0.408 factor: it converts radiation from MJ m−2 into millimetres of evaporated water; without it the radiation term is 2.45 times too large.

06 lépés step

Kiválasztom a megfelelő nap időjárását I pick the weather of the right day

Ezt az eredeti TIF-ek reprodukálásakor vettem észre. A script így keresi az ERA5-képet: filterDate(date, date.advance(1,'day')).first(), ahol a date a műhold áthaladásának időpontja (~09:50 UTC). A napi ERA5-kép viszont 00:00-kor kezdődik, így a szűrőablakból kiesik – a first() a következő nap képét adja.

Váltogasd az ábrán a két módot. Az „Eredeti script” mód szándékosan megtartja ezt a viselkedést (így pixelre egyezik a régi eredményekkel), a „Javított FAO-56” mód a felvétel napjának időjárásával számol.

I noticed this while reproducing the original TIFs. The script looks up ERA5 like this: filterDate(date, date.advance(1,'day')).first(), where date is the time of the satellite overpass (~09:50 UTC). But the daily ERA5 image starts at 00:00, so it falls outside the window – first() returns the image of the next day.

Switch between the two modes in the figure. The “Original script” mode keeps this behaviour on purpose (so it matches the old results pixel by pixel); the “Corrected FAO-56” mode uses the weather of the acquisition day.

07 lépés step

Összeszorzom: ETa = Kc × ET0 I multiply: ETa = Kc × ET0

A végeredmény pixelenként a Kc és az ET0 szorzata, mm/nap-ban. A Kc 10 m-es (Sentinel-2), az ET0 viszont ~9 km-es (ERA5-Land) – egy ERA5-cellára ~800 × 800 Sentinel-pixel jut.

Az eredeti script legközelebbi-szomszéd mintavételt használt: az ERA5-cellák határán a térképen éles lépcsők jelennek meg. A javított módban bilineárisan interpolálom az ERA5-mezőt, így az időjárás simán változik a területen. Kapcsold át az ábrát!

The result is, pixel by pixel, the product of Kc and ET0 in mm/day. Kc is 10 m (Sentinel-2) but ET0 is ~9 km (ERA5-Land) – one ERA5 cell holds ~800 × 800 Sentinel pixels.

The original script used nearest-neighbour sampling: sharp steps appear on the map along the ERA5 cell borders. In the corrected mode I interpolate the ERA5 field bilinearly, so the weather changes smoothly across the area. Toggle the figure!

08 lépés step

Statisztikát számolok és GeoTIFF-et adok I compute statistics and hand over GeoTIFFs

Ha poligont rajzolsz, az Earth Engine-ben pontosan a GeoTIFF-ek pixelrácsán számolom ki az ETa átlagát, mediánját, szórását, minimumát, maximumát és a pixelszámot – minden napra. Egy pixel akkor számít, ha a középpontja a poligonon belül van (ábra). A medián 0,001 mm/nap pontosságú.

A letöltésnél a nagy területet legfeljebb 8 megapixeles csempékben kérem le az Earth Engine-ből, majd egyetlen tömörített, float32 GeoTIFF-be fűzöm (NoData = −9999). A ZIP-be egy ETa_run_info.json is kerül, minden beállítással, hogy a futás bármikor megismételhető legyen.

When you draw a polygon, I compute the mean, median, standard deviation, minimum, maximum and pixel count of ETa in Earth Engine on exactly the pixel grid of the GeoTIFFs – for every date. A pixel counts when its centre is inside the polygon (figure). The median is exact to 0.001 mm/day.

For downloads I request large areas from Earth Engine in tiles of at most 8 megapixels and merge them into one compressed float32 GeoTIFF (NoData = −9999). The ZIP also gets an ETa_run_info.json with every setting, so the run can be repeated at any time.

Javítások Corrections

Hol és miért javítottam az eredeti scripten? Where and why I corrected the original script

Mindegyik javítás a „Javított FAO-56” módban él. Az „Eredeti script” mód szándékosan mindet kihagyja, hogy a régi eredmények pontosan reprodukálhatók legyenek. Every correction is active in the “Corrected FAO-56” mode. The “Original script” mode leaves them all out on purpose, so the old results can be reproduced exactly.

1

Hiányzó 0,408-as szorzó Missing 0.408 factor

GEE_code.txt · 5. calcET0
eredeti original
ET0 = delta.multiply(Rn)
  .add(...)
javított corrected
ET0 = 0.408 · Δ · Rn
  + γ·900/(T+273)·u2·(es−ea)

A FAO-56 egyenletben a 0,408 alakítja át a sugárzási energiát (MJ m−2 nap−1) párolgó vízmagassággá (mm nap−1) – ez a víz párolgáshőjének reciproka. Nélküle a sugárzási tag 1/0,408 = 2,45-ször nagyobb. In the FAO-56 equation 0.408 converts radiation energy (MJ m−2 day−1) into evaporated water depth (mm day−1) – it is the reciprocal of the latent heat of vaporisation. Without it the radiation term is 1/0.408 = 2.45 times larger.

A legnagyobb hatású hiba: önmagában ~2× túlbecslés. The largest error: about 2× overestimate on its own.
2

Nettó sugárzás hosszúhullám nélkül Net radiation without longwave

GEE_code.txt · 5. Rn
eredeti original
Rn = surface_net_solar_radiation_sum
javított corrected
Rn = 0.77 · Rs
   + surface_net_thermal_radiation_sum

A script csak a rövidhullámú (napsugárzás) mérleget vette, a felszín hosszúhullámú kisugárzását (éjjel-nappal napi 3–6 MJ m−2) nem vonta le. A FAO-56 referencia-gyep albedója 0,23, ezért a rövidhullámú tagot 0,77·Rs-ből számolom. The script only took the shortwave (solar) balance and did not subtract the longwave emission of the surface (3–6 MJ m−2 per day). The FAO-56 reference grass has albedo 0.23, so I compute the shortwave term as 0.77·Rs.

Rn kb. 20–30%-kal kisebb. Rn about 20–30% smaller.
3

10 m-es szél 2 m-es helyett 10 m wind used as 2 m wind

GEE_code.txt · 5. u2
eredeti original
u2 = sqrt(u10² + v10²)   // napi átlag komponensek
javított corrected
u2 = mean(|u10,v10| óránként) · 4.87 / ln(67.8·10 − 5.42)

A FAO-56 a 2 m-en mért szélsebességet várja; 10 m-en a szél ~25%-kal erősebb. Ráadásul a napi átlagos u- és v-komponensből számolt sebesség kioltja az irányváltó szelet. Az óránkénti sebességek átlagát veszem, és a FAO-56 logaritmikus profiljával átszámolom 2 m-re (×0,748). FAO-56 expects wind speed at 2 m; at 10 m the wind is ~25% stronger. Also, a speed computed from the daily mean u and v components cancels out changing wind directions. I take the mean of the hourly speeds and convert it to 2 m with the logarithmic profile of FAO-56 (×0.748).

Kisebb, de rendszeres hatás a légköri tagon. Smaller but systematic effect on the aerodynamic term.
4

Telítési gőznyomás és pszichrometrikus állandó Saturation vapour pressure and psychrometric constant

GEE_code.txt · 5. es, γ
eredeti original
es = e°(Tátlag)
gamma = 0.066
javított corrected
es = (e°(Tmax) + e°(Tmin)) / 2
γ = 0.000665 · P

Az e°(T) görbe nemlineáris, ezért a napi átlaghőmérsékletből számolt es alulbecsli a telítési hiányt; a FAO-56 a Tmax és Tmin értékekből számolt e° átlagát írja elő. A γ a légnyomással arányos – az ERA5-Land felszíni nyomásából számolom. The e°(T) curve is non-linear, so es from the mean temperature underestimates the vapour pressure deficit; FAO-56 prescribes the mean of e° at Tmax and Tmin. γ is proportional to air pressure – I compute it from the ERA5-Land surface pressure.

Nyáron +5–15% a légköri tagban, a γ-hatás kicsi. In summer +5–15% in the aerodynamic term; the γ effect is small.
5

A következő nap időjárása The weather of the next day

GEE_code.txt · 4. getEra5Image
eredeti original
filterDate(date, date.advance(1,'day'))
// date = áthaladás ~09:50 UTC
javított corrected
filterDate(nap 00:00, nap+1 00:00)

Az áthaladás időpontjától induló egynapos ablak kihagyja a felvétel napjának ERA5-képét, és a következő napét adja vissza. A javított mód a felvétel napjával számol. The one-day window starting at the overpass time skips the ERA5 image of the acquisition day and returns the next day's. The corrected mode uses the acquisition day.

Napról napra ±3–15%; az eredeti TIF-ek reprodukálásakor derült ki. ±3–15% day to day; found while reproducing the original TIFs.
6

Sentinel-2 gyűjtemény, felhőmaszk, mozaik Sentinel-2 collection, cloud mask, mosaic

GEE_code.txt · 2., 7.
eredeti original
COPERNICUS/S2_SR  (elavult)
csak CLOUDY_PIXEL_PERCENTAGE
az első 10 felvétel
javított corrected
COPERNICUS/S2_SR_HARMONIZED
+ SCL pixel-felhőmaszk
minden nap egy mozaik

Az elavult gyűjteményben 2022 óta +1000-es eltolás torzítja az NDVI-t; a felhős pixelek maszk nélkül is értéket kaptak; és a script csak az első 10 felvételt exportálta, csempénként külön. In the deprecated collection a +1000 offset has biased NDVI since 2022; cloudy pixels got values without a mask; and the script exported only the first 10 scenes, one tile at a time.

Helyes NDVI → helyes Kc; a felhős pixelek üresek. Correct NDVI → correct Kc; cloudy pixels are empty.
+

Lépcsős ERA5 helyett sima mező Smooth ERA5 field instead of steps

GEE_code.txt · 5. (mintavétel)
eredeti original
ERA5 legközelebbi szomszéd
(~9 km-es blokkok)
javított corrected
ERA5 .resample('bilinear')

Ez nem hiba a szó szoros értelmében, de a 9 km-es cellahatárokon mesterséges ugrások jelentek meg az ETa-térképen, amelyek a poligon-statisztikákat is torzíthatják. Not an error in the strict sense, but artificial jumps appeared on the ETa map at the 9 km cell borders, which can also bias polygon statistics.

Nincs cellahatár-lépcső a térképen. No cell-border steps on the map.
Ellenőrzés Validation

Az „Eredeti script” mód pixelre reprodukálja a régi eredményeket The “Original script” mode reproduces the old results pixel by pixel

Két korábbi, a scripttel készült GeoTIFF-et (2023-04-22 és 2023-07-16) vetettem össze az app kimenetével: a 14,9 millió közös pixel mindegyike azonos, a legnagyobb eltérés 0. A javított mód ugyanezekre a napokra 2,5–3-szor kisebb, fizikailag reális értéket ad. I compared two earlier GeoTIFFs made with the script (2023-04-22 and 2023-07-16) with the app's output: all 14.9 million common pixels are identical, the largest difference is 0. The corrected mode gives 2.5–3 times smaller, physically realistic values for the same days.

Források References

  • Allen, R.G., Pereira, L.S., Raes, D., Smith, M. (1998): Crop evapotranspiration – Guidelines for computing crop water requirements. FAO Irrigation and Drainage Paper 56, Rome. (eq. 6, 8, 11–14, 38–39, 47)
  • Muñoz-Sabater, J. et al. (2021): ERA5-Land: a state-of-the-art global reanalysis dataset for land applications. Earth Syst. Sci. Data 13, 4349–4383. · Earth Engine: ECMWF/ERA5_LAND/DAILY_AGGR, ECMWF/ERA5_LAND/HOURLY
  • Copernicus Sentinel-2 L2A, Earth Engine: COPERNICUS/S2_SR_HARMONIZED (processing baseline 04.00 offset since 2022-01-25)
ETa Mapper v1.0 · Sentinel-2 (Copernicus) · ERA5-Land (ECMWF / C3S) · Google Earth Engine · HUN-REN CSFK FGI
file