|
| 1 | +#!/usr/bin/env python3 |
| 2 | +# -*- coding: utf-8 -*- |
| 3 | +""" |
| 4 | +Created on Tue Jul 22 10:51:59 2025 |
| 5 | +
|
| 6 | +@author: Hisham Eldardiry |
| 7 | +""" |
| 8 | + |
| 9 | +import ee |
| 10 | +import pandas as pd |
| 11 | +from datetime import datetime, timedelta |
| 12 | + |
| 13 | +# Authenticate and initialize Earth Engine |
| 14 | +ee.Initialize() |
| 15 | + |
| 16 | + |
| 17 | +# --- CONFIGURATION --- |
| 18 | +lat, lon = 11.9189, 107.9069 |
| 19 | +start_date = '2022-01-01' |
| 20 | +end_date = '2022-12-31' |
| 21 | + |
| 22 | +point = ee.Geometry.Point([lon, lat]) |
| 23 | + |
| 24 | +# --- Generate daily list of dates --- |
| 25 | +start_dt = datetime.strptime(start_date, "%Y-%m-%d") |
| 26 | +end_dt = datetime.strptime(end_date, "%Y-%m-%d") |
| 27 | +date_list = [(start_dt + timedelta(days=i)).strftime("%Y-%m-%d") |
| 28 | + for i in range((end_dt - start_dt).days + 1)] |
| 29 | + |
| 30 | + |
| 31 | +# --- CHIRPS Daily Precipitation --- |
| 32 | +chirps = ee.ImageCollection("UCSB-CHG/CHIRPS/DAILY") \ |
| 33 | + .filterDate(start_date, end_date) \ |
| 34 | + .filterBounds(point) |
| 35 | + |
| 36 | +# --- ERA5-Land Hourly --- |
| 37 | +era5 = ee.ImageCollection("ECMWF/ERA5_LAND/HOURLY") \ |
| 38 | + .filterDate(start_date, end_date) \ |
| 39 | + .filterBounds(point) \ |
| 40 | + .select(['temperature_2m', 'u_component_of_wind_10m', 'v_component_of_wind_10m']) |
| 41 | + |
| 42 | +# %% Add derived bands: temperature in C, wind speed |
| 43 | +def process_hourly(img): |
| 44 | + temp_c = img.select('temperature_2m').subtract(273.15).rename('temp_C') |
| 45 | + u = img.select('u_component_of_wind_10m') |
| 46 | + v = img.select('v_component_of_wind_10m') |
| 47 | + wind = u.pow(2).add(v.pow(2)).sqrt().rename('wind_speed') |
| 48 | + return img.addBands([temp_c, wind]) |
| 49 | + |
| 50 | +era5 = era5.map(process_hourly) |
| 51 | + |
| 52 | + |
| 53 | + |
| 54 | + |
| 55 | +# %% Function to create daily combined image |
| 56 | +def daily_reduce(date_str): |
| 57 | + date = ee.Date(date_str) |
| 58 | + next_date = date.advance(1, 'day') |
| 59 | + |
| 60 | + # Handle missing CHIRPS |
| 61 | + chirps_filtered = chirps.filterDate(date, next_date) |
| 62 | + chirps_img = ee.Image(0).rename('precipitation') # default image |
| 63 | + chirps_img = ee.Image(ee.Algorithms.If( |
| 64 | + chirps_filtered.size().gt(0), |
| 65 | + chirps_filtered.mean().select('precipitation'), |
| 66 | + chirps_img |
| 67 | + )) |
| 68 | + |
| 69 | + # ERA5: tmin, tmax, wind |
| 70 | + daily_imgs = era5.filterDate(date, next_date) |
| 71 | + tmin = daily_imgs.select('temp_C').min() |
| 72 | + tmax = daily_imgs.select('temp_C').max() |
| 73 | + wind = daily_imgs.select('wind_speed').mean() |
| 74 | + |
| 75 | + combined = chirps_img.addBands([tmax.rename('tmax'), tmin.rename('tmin'), wind.rename('wind')]) |
| 76 | + return combined.set('date', date.format('YYYY-MM-dd')) |
| 77 | + |
| 78 | +# Create image collection for daily combined records |
| 79 | +daily_images = ee.ImageCollection([daily_reduce(d) for d in date_list]) |
| 80 | + |
| 81 | +# %% Extract data as list of features |
| 82 | +def extract_record(img): |
| 83 | + values = img.reduceRegion( |
| 84 | + reducer=ee.Reducer.first(), |
| 85 | + geometry=point, |
| 86 | + scale=5000 |
| 87 | + ) |
| 88 | + return ee.Feature(None, values).set('date', img.get('date')) |
| 89 | + |
| 90 | +feature_collection = ee.FeatureCollection(daily_images.map(extract_record)) |
| 91 | +features = feature_collection.getInfo()['features'] |
| 92 | + |
| 93 | + |
| 94 | +# %% Convert to DataFrame |
| 95 | +df = pd.DataFrame([ |
| 96 | + { |
| 97 | + 'date': f['properties']['date'], |
| 98 | + 'precip': f['properties'].get('precipitation'), |
| 99 | + 'tmax': f['properties'].get('tmax'), |
| 100 | + 'tmin': f['properties'].get('tmin'), |
| 101 | + 'wind': f['properties'].get('wind') |
| 102 | + } |
| 103 | + for f in features |
| 104 | +]) |
| 105 | + |
| 106 | +# %% Save Output |
| 107 | +df[['precip', 'tmax', 'tmin', 'wind']].to_csv( |
| 108 | + f"data_{lat:.4f}_{lon:.4f}", sep=' ', index=False, header=False |
| 109 | +) |
| 110 | + |
| 111 | +print(f"Saved forcing file for lat {lat}, lon {lon}") |
0 commit comments