— This UDF implements the Gumbel distribution formula to estimate the event magnitude — for a given return period based on the sample mean (xbar) and standard deviation (sigma). CREATE TEMP FUNCTION CalculateReturnPeriod(period INT64, xbar FLOAT64, sigma FLOAT64) RETURNS FLOAT64 AS ( ROUND(-LOG(-LOG(1 – (1 / period))) * sigma * .7797 + xbar – (.45 * sigma), 4) ); WITH — Step 1: Define the analysis areas. — This CTE selects a specific subset of 10 major German cities. — ST_SIMPLIFY is used to reduce polygon complexity, improving query performance. Counties AS ( FROM bigquery-public-data.overture_maps.division_area |> WHERE country = ‘DE’ AND subtype = ‘county’ AND names.primary IN ( ‘München’, ‘Köln’, ‘Frankfurt am Main’, ‘Stuttgart’, ‘Düsseldorf’) |> SELECT id, names.primary AS county, ST_SIMPLIFY(geometry,500) AS geometry ), — Step 2: Define the time periods for comparison. — These two 30-year, overlapping epochs will be used to assess recent changes. Epochs AS ( FROM UNNEST([ STRUCT(‘e1’ AS epoch, 1980 AS start_year, 2010 AS end_year), STRUCT(‘e2’ AS epoch, 1994 AS start_year, 2024 AS end_year)]) ), — Step 3: Define the return periods to calculate. ReturnPeriods AS ( FROM UNNEST([10,25,50,100,500]) AS years |> SELECT * ), — Step 4: Select the relevant image data from the Earth Engine catalog. — Replace YOUR_CLOUD_PROJECT with your relevant Cloud Project ID. Images AS ( FROM YOUR_CLOUD_PROJECT.era5_land_daily_aggregated.climate |> WHERE year BETWEEN 1980 AND 2024 |> SELECT id AS img_id, start_datetime AS img_date ) — Step 5: Begin the main data processing pipeline. — This creates a processing task for every combination of a day and a county. FROM Images |> CROSS JOIN Counties — Step 6: Perform zonal statistics using Earth Engine. — ST_REGIONSTATS calculates the mean precipitation for each county for each day. |> SELECT img_id, Counties.id AS county_id, EXTRACT(YEAR FROM img_date) AS year, ST_REGIONSTATS(geometry, img_id, ‘total_precipitation_sum’) AS areal_precipitation_stats — Step 7: Find the annual maximum precipitation. — This aggregates the daily results to find the single wettest day for each county within each year. |> AGGREGATE SUM(areal_precipitation_stats.count) AS pixels_examined, MAX(areal_precipitation_stats.mean) AS yearly_max_1day_precip, ANY_VALUE(areal_precipitation_stats.area) AS pixel_area GROUP BY county_id, year — Step 8: Calculate statistical parameters for each epoch. — Joins the annual maxima to the epoch definitions and then calculates the — average and standard deviation required for the Gumbel distribution formula. |> JOIN Epochs ON year BETWEEN start_year AND end_year |> AGGREGATE AVG(yearly_max_1day_precip * 1e3) AS avg_yearly_max_1day_precip, STDDEV(yearly_max_1day_precip * 1e3) AS stddev_yearly_max_1day_precip, GROUP BY county_id, epoch — Step 9: Calculate the return period precipitation values. — Applies the UDF to the calculated statistics for each return period. — This assumes they are the same function. |> CROSS JOIN ReturnPeriods rp |> EXTEND CalculateReturnPeriod(rp.years, avg_yearly_max_1day_precip, stddev_yearly_max_1day_precip) AS est_max_1day_precip |> DROP avg_yearly_max_1day_precip, stddev_yearly_max_1day_precip — Step 10: Pivot the data to create columns for each epoch and return period. — The first PIVOT transforms rows for ‘e1’ and ‘e2’ into columns for direct comparison. |> PIVOT (ANY_VALUE(est_max_1day_precip) AS est_max_1day_precip FOR epoch IN (‘e1’, ‘e2’)) — The second PIVOT transforms rows for each return period (10, 25, etc.) into columns. |> PIVOT (ANY_VALUE(est_max_1day_precip_e1) AS e1, ANY_VALUE(est_max_1day_precip_e2) AS e2 FOR years IN (10, 25, 50, 100, 500)) — Step 11: Re-attach county names and geometries for the final output. |> LEFT JOIN Counties ON county_id = Counties.id — Step 12: Calculate the final difference between the two epochs. — This creates the delta values that show the change in precipitation magnitude for each return period. |> SELECT county, Counties.geometry, e2_10 – e1_10 AS est_10yr_max_1day_precip_delta, e2_25 – e1_25 AS est_25yr_max_1day_precip_delta, e2_50 – e1_50 AS est_50yr_max_1day_precip_delta, e2_100 – e1_100 AS est_100yr_max_1day_precip_delta, e2_500 – e1_500 AS est_500yr_max_1day_precip_delta |> ORDER BY county
