-
Notifications
You must be signed in to change notification settings - Fork 93
Description
Hi everyone,
I have got an error that I don't arrive to correct despite all my efforts.
Here is my code : the final aim is to compare the yield for an AgriPV and a not AgriPV field. So, first I am collecting wheater date from NASA POWER. Then, I convert all datas in the dataframe requiered by aquacrop. But finally my model do not run and I can't retrieve yield data. Can you help me please ? Thank you !
Here is my code, the error message is about run_simulation() fonction and the success = model.run_model(till_termination=True) command
`import pandas as pd
import numpy as np
import requests
from math import radians, tan, sin, cos
from tqdm import tqdm
from pvlib.solarposition import get_solarposition
import matplotlib.pyplot as plt
import requests
import pandas as pd
==================== 1. Weather and hourly solar radiation data ====================
def get_nasa_power_data_hourly(lat, lon, start='2021-01-01', end='2021-12-31'):
url = "https://power.larc.nasa.gov/api/temporal/hourly/point"
params = {
"parameters": "ALLSKY_SFC_SW_DWN,ALLSKY_SFC_SW_DIFF,ALLSKY_SFC_SW_DNI,ALLSKY_SRF_ALB",
"community": "RE",
"longitude": lon,
"latitude": lat,
"start": start.replace("-", ""),
"end": end.replace("-", ""),
"format": "JSON",
"time-standard": "UTC"
}
print("📡 Requesting NASA POWER data...")
response = requests.get(url, params=params)
if response.status_code != 200:
print("❌ NASA POWER request error:", response.status_code)
print("↪ Raw response:", response.text)
raise RuntimeError("NASA POWER API error")
data = response.json()
param = data.get("properties", {}).get("parameter", {})
ghi = param.get("ALLSKY_SFC_SW_DWN", {})
dhi = param.get("ALLSKY_SFC_SW_DIFF", {})
dni = param.get("ALLSKY_SFC_SW_DNI", {})
albedo = param.get("ALLSKY_SRF_ALB", {})
if not (ghi and dhi and dni and albedo):
print("⚠️ Missing data for some parameters.")
return None
# Extract common timestamps
timestamps = sorted(set(ghi.keys()) & set(dhi.keys()) & set(dni.keys()) & set(albedo.keys()))
df = pd.DataFrame({
"date": pd.to_datetime(timestamps, format="%Y%m%d%H", utc=True),
"GHI": [ghi[t] for t in timestamps],
"DHI": [dhi[t] for t in timestamps],
"DNI": [dni[t] for t in timestamps],
"ALB": [albedo[t] for t in timestamps]
})
df["GHI"] = df["GHI"].astype(float)
df["DHI"] = df["DHI"].astype(float)
df["DNI"] = df["DNI"].astype(float)
df["ALB"] = df["ALB"].astype(float)
df["day_of_year"] = df["date"].dt.dayofyear
print(f"✅ {len(df)} hourly values loaded.")
return df
==================== 2. Solar geometry ====================
def compute_solar_geometry_hourly(df, lat, lon):
times = df.index.tz_convert(None)
solpos = get_solarposition(times, lat, lon)
df['solar_zenith'] = solpos['apparent_zenith'].values
df['solar_azimuth'] = solpos['azimuth'].values
df['solar_elevation'] = 90 - df['solar_zenith']
return df
==================== 3. Shadow simulation ====================
def simulate_shadow_map(df, panel_tilt, row_spacing, panel_height, leg_height, grid_size=10):
dx = 1.0 # m/pixel
n_pixels = grid_size
field_size = (n_pixels, n_pixels)
irradiance_maps = np.zeros((len(df), *field_size))
row_interval = int(row_spacing / dx)
panel_rows = list(range(0, n_pixels, row_interval))
for t, (index, row) in tqdm(enumerate(df.iterrows()), total=len(df), desc="Shadow calculation"):
elevation = row["solar_elevation"]
azimuth = row["solar_azimuth"]
if elevation <= 0:
continue
tilt_rad = radians(panel_tilt)
h_panel = leg_height + panel_height * sin(tilt_rad)
shadow_length = h_panel / tan(radians(elevation))
shadow_pixels = int(round(shadow_length / dx))
azimuth_rad = radians(azimuth)
dx_shadow = -sin(azimuth_rad)
dy_shadow = -cos(azimuth_rad)
shadow_map = np.ones(field_size)
for y_row in panel_rows:
for x in range(n_pixels):
for s in range(1, shadow_pixels + 1):
xx = int(round(x + s * dx_shadow))
yy = int(round(y_row + s * dy_shadow))
if 0 <= xx < n_pixels and 0 <= yy < n_pixels:
shadow_map[xx, yy] *= 0.2
irradiance_maps[t] = shadow_map * row["GHI"]
return irradiance_maps
==================== 4. Daily aggregation ====================
def aggregate_irradiance(irradiance_maps, df):
hourly_means = np.mean(irradiance_maps.reshape(len(df), -1), axis=1)
df = df.copy()
df['irradiance_corrected'] = hourly_means
df['date_only'] = df.index.date
df_daily = df.groupby('date_only')['irradiance_corrected'].sum().reset_index()
df_daily.columns = ['date', 'SolarRadiation']
return df_daily
==================== 5. Daily weather data ====================
def get_nasa_weather_daily(lat, lon, start, end):
url = "https://power.larc.nasa.gov/api/temporal/daily/point"
params = {
"parameters": "T2M_MIN,T2M_MAX,PRECTOTCORR",
"community": "AG",
"longitude": lon,
"latitude": lat,
"start": start.replace("-", ""),
"end": end.replace("-", ""),
"format": "JSON"
}
response = requests.get(url, params=params)
data = response.json()["properties"]["parameter"]
dates = pd.date_range(start=start, end=end)
df = pd.DataFrame({
"date": dates,
"MinTemp": [data["T2M_MIN"][d.strftime("%Y%m%d")] for d in dates],
"MaxTemp": [data["T2M_MAX"][d.strftime("%Y%m%d")] for d in dates],
"Precipitation": [data["PRECTOTCORR"][d.strftime("%Y%m%d")] for d in dates],
})
return df
==================== 6. Merge weather and radiation data ====================
def prepare_aquacrop_input(lat, lon, start, end, shadow=True):
df_hourly = get_nasa_power_data_hourly(lat, lon, start, end)
df_hourly.set_index("date", inplace=True)
df_hourly = compute_solar_geometry_hourly(df_hourly, lat, lon)
if shadow:
irradiance_maps = simulate_shadow_map(df_hourly, panel_tilt=25, row_spacing=5, panel_height=1, leg_height=1)
irradiance_daily = aggregate_irradiance(irradiance_maps, df_hourly)
else:
df_hourly['date_only'] = df_hourly.index.date
irradiance_daily = df_hourly.groupby('date_only')['GHI'].sum().reset_index()
irradiance_daily.columns = ['date', 'SolarRadiation']
meteo = get_nasa_weather_daily(lat, lon, start, end)
meteo['date'] = pd.to_datetime(meteo['date'])
irradiance_daily['date'] = pd.to_datetime(irradiance_daily['date'])
meteo = meteo.merge(irradiance_daily, on='date')
meteo['ET0'] = 0.0023 * (meteo["MaxTemp"] - meteo["MinTemp"])**0.5 * (meteo["MaxTemp"] + meteo["MinTemp"]) / 2 + 0.408 * meteo["SolarRadiation"] / 2.45
# 🔧 Rename columns for AquaCrop
meteo = meteo.rename(columns={
'date': 'Date',
'ET0': 'ReferenceET'})
meteo = meteo[['Date', 'MinTemp', 'MaxTemp', 'Precipitation', 'ReferenceET']]
meteo["MinTemp"] = pd.to_numeric(meteo["MinTemp"], errors="coerce")
meteo["MaxTemp"] = pd.to_numeric(meteo["MaxTemp"], errors="coerce")
meteo["Precipitation"] = pd.to_numeric(meteo["Precipitation"], errors="coerce")
meteo["ReferenceET"] = pd.to_numeric(meteo["ReferenceET"], errors="coerce")
print("🧪 Final data ready for AquaCrop:", meteo.columns)
return meteo
==================== 7. Export to AquaCrop format ====================
def export_to_aquacrop_format(df, path):
df_out = df[['Date', 'MinTemp', 'MaxTemp', 'Precipitation', 'ReferenceET']].copy()
df_out.columns = ['Date', 'Tmin', 'Tmax', 'P', 'ET0']
df_out['Date'] = df_out['Date'].dt.strftime('%Y-%m-%d')
df_out.to_csv(path, index=False)
==================== 8. AquaCrop simulation using pyAquaCrop ====================
from aquacrop import AquaCropModel, Soil, Crop, InitialWaterContent
def run_simulation(meteo_df, name):
# meteo_df must contain these columns: Date, MinTemp, MaxTemp, Precipitation, ReferenceET
print("\n🛠️ Checking data types for MinTemp / MaxTemp:")
print(meteo_df[["MinTemp", "MaxTemp"]].applymap(type).value_counts())
print("\n🔎 Looking for Timestamp in MinTemp:")
print(meteo_df[meteo_df["MinTemp"].apply(lambda x: isinstance(x, pd.Timestamp))])
for col in ['MinTemp', 'MaxTemp']:
wrong_types = meteo_df[~meteo_df[col].apply(lambda x: isinstance(x, (int, float)))]
if not wrong_types.empty:
print(f"\n❌ Problem detected in {col}:")
print(wrong_types.head())
raise TypeError(f"Column {col} contains non-numeric types (e.g., {type(wrong_types.iloc[0][col])})")
print("\n📋 Unique types per column:")
for col in ['MinTemp', 'MaxTemp', 'Precipitation', 'ReferenceET']:
meteo_df[col] = pd.to_numeric(meteo_df[col], errors="coerce")
print(f"{col}: {meteo_df[col].apply(type).value_counts()}")
start_date = pd.to_datetime(meteo_df["Date"]).iloc[0].strftime("%Y/%m/%d")
end_date = pd.to_datetime(meteo_df["Date"]).iloc[-1].strftime("%Y/%m/%d")
print(f"Start: {start_date}, End: {end_date}, Type: {type(start_date)}")
soil = Soil("SandyLoam")
crop = Crop("Maize", planting_date="04/01")
init_wc = InitialWaterContent("FC")
model = AquaCropModel(
start_date,
end_date,
meteo_df,
soil,
crop,
init_wc
)
success = model.run_model(till_termination=True)
if not success:
print("Simulation failed")
return None
return model
def validate_and_clean_meteo_df(df):
for col in ['MinTemp', 'MaxTemp', 'Precipitation', 'ReferenceET']:
df[col] = pd.to_numeric(df[col], errors='coerce')
for col in ['MinTemp', 'MaxTemp', 'Precipitation', 'ReferenceET']:
non_numeric = df[~df[col].apply(lambda x: isinstance(x, (int, float)))]
if not non_numeric.empty:
print(f"Non-numeric values found in {col}:")
print(non_numeric)
raise ValueError(f"Non-numeric values found in {col}")
df = df.dropna(subset=['MinTemp', 'MaxTemp', 'Precipitation', 'ReferenceET'])
return df
==================== 9. usage ====================
lat, lon = 45.0, 3.0
start = "2021-01-01"
end = "2021-12-31"
meteo_shadow = prepare_aquacrop_input(lat, lon, start, end, shadow=True)
meteo_clear = prepare_aquacrop_input(lat, lon, start, end, shadow=False)
print(meteo_shadow.dtypes)
print(meteo_shadow.head(10))
export_to_aquacrop_format(meteo_shadow, "meteo_with_panels.csv")
export_to_aquacrop_format(meteo_clear, "meteo_without_panels.csv")
print("Types in meteo_shadow['MinTemp']:")
print(meteo_shadow["MinTemp"].apply(type).value_counts())
print("Types in meteo_shadow['MaxTemp']:")
print(meteo_shadow["MaxTemp"].apply(type).value_counts())
bad_rows = meteo_shadow[meteo_shadow["MinTemp"].apply(lambda x: isinstance(x, pd.Timestamp))]
print("Rows with Timestamp in MinTemp:")
print(bad_rows)
meteo_shadow["MinTemp"] = pd.to_numeric(meteo_shadow["MinTemp"], errors="coerce")
meteo_shadow["MaxTemp"] = pd.to_numeric(meteo_shadow["MaxTemp"], errors="coerce")
meteo_shadow = meteo_shadow.dropna(subset=["MinTemp", "MaxTemp"])
meteo_shadow = validate_and_clean_meteo_df(meteo_shadow)
meteo_clear = validate_and_clean_meteo_df(meteo_clear)
model_shadow = run_simulation(meteo_shadow, "With panels")
model_clear = run_simulation(meteo_clear, "Without panels")
`
Here is the exact error message
`📡 Requesting NASA POWER data...
✅ 8760 hourly values loaded.
Shadow calculation: 100%|██| 8760/8760 [00:03<00:00, 2822.81it/s]
🧪 Final data ready for AquaCrop: Index(['Date', 'MinTemp', 'MaxTemp', 'Precipitation', 'ReferenceET'], dtype='object')
📡 Requesting NASA POWER data...
✅ 8760 hourly values loaded.
🧪 Final data ready for AquaCrop: Index(['Date', 'MinTemp', 'MaxTemp', 'Precipitation', 'ReferenceET'], dtype='object')
Date datetime64[ns]
MinTemp float64
MaxTemp float64
Precipitation float64
ReferenceET float64
dtype: object
Date MinTemp MaxTemp Precipitation ReferenceET
0 2021-01-01 -10.41 -1.60 1.15 67.094424
1 2021-01-02 -7.29 -3.17 0.18 43.354624
2 2021-01-03 -9.72 -3.03 0.44 78.721383
3 2021-01-04 -9.56 -1.27 3.44 79.152975
4 2021-01-05 -10.59 -1.86 0.43 50.847191
5 2021-01-06 -10.81 -2.99 0.10 48.384140
6 2021-01-07 -10.07 -2.62 0.22 75.721669
7 2021-01-08 -9.55 -2.73 0.09 98.918739
8 2021-01-09 -7.45 -3.58 0.03 103.869624
9 2021-01-10 -9.10 -4.25 0.03 97.965490
Types in meteo_shadow['MinTemp']:
MinTemp
<class 'float'> 365
Name: count, dtype: int64
Types in meteo_shadow['MaxTemp']:
MaxTemp
<class 'float'> 365
Name: count, dtype: int64
Rows with Timestamp in MinTemp:
Empty DataFrame
Columns: [Date, MinTemp, MaxTemp, Precipitation, ReferenceET]
Index: []
🛠️ Checking data types for MinTemp / MaxTemp:
c:\Users\omare\OneDrive\Documents\UNSW\Research\RECHERCHE\VScode\projet_stics\mon_script_aquacrop.py:185: FutureWarning: DataFrame.applymap has been deprecated. Use DataFrame.map instead.
print(meteo_df[["MinTemp", "MaxTemp"]].applymap(type).value_counts())
MinTemp MaxTemp
<class 'float'> <class 'float'> 365
Name: count, dtype: int64
🔎 Looking for Timestamp in MinTemp:
Empty DataFrame
Columns: [Date, MinTemp, MaxTemp, Precipitation, ReferenceET]
Index: []
📋 Unique types per column:
MinTemp: MinTemp
<class 'float'> 365
Name: count, dtype: int64
MaxTemp: MaxTemp
<class 'float'> 365
Name: count, dtype: int64
Precipitation: Precipitation
<class 'float'> 365
Name: count, dtype: int64
ReferenceET: ReferenceET
<class 'float'> 365
Name: count, dtype: int64
Start: 2021/01/01, End: 2021/12/31, Type: <class 'str'>
Traceback (most recent call last):
File "c:\Users\omare\OneDrive\Documents\UNSW\Research\RECHERCHE\VScode\projet_stics\mon_script_aquacrop.py", line 268, in
model_shadow = run_simulation(meteo_shadow, "With panels")
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "c:\Users\omare\OneDrive\Documents\UNSW\Research\RECHERCHE\VScode\projet_stics\mon_script_aquacrop.py", line 219, in run_simulation
success = model.run_model(till_termination=True)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "C:\Users\omare\OneDrive\Documents\UNSW\Research\RECHERCHE\VScode.venv\Lib\site-packages\aquacrop\core.py", line 279, in run_model
) = self._perform_timestep()
^^^^^^^^^^^^^^^^^^^^^^^^
File "C:\Users\omare\OneDrive\Documents\UNSW\Research\RECHERCHE\VScode.venv\Lib\site-packages\aquacrop\core.py", line 325, in _perform_timestep
new_cond, param_struct, outputs = solution_single_time_step(
^^^^^^^^^^^^^^^^^^^^^^^^^^
File "C:\Users\omare\OneDrive\Documents\UNSW\Research\RECHERCHE\VScode.venv\Lib\site-packages\aquacrop\timestep\run_single_timestep.py", line 136, in solution_single_time_step
gdd = growing_degree_day(
^^^^^^^^^^^^^^^^^^^
File "C:\Users\omare\OneDrive\Documents\UNSW\Research\RECHERCHE\VScode.venv\Lib\site-packages\aquacrop\solution\growing_degree_day.py", line 58, in growing_degree_day
temp_min = min(temp_min, Tupp)
^^^^^^^^^^^^^^^^^^^
TypeError: '<' not supported between instances of 'float' and 'Timestamp'
`