diff --git a/quartz_solar_forecast/eval/pv.py b/quartz_solar_forecast/eval/pv.py index 091734824..529596bef 100644 --- a/quartz_solar_forecast/eval/pv.py +++ b/quartz_solar_forecast/eval/pv.py @@ -2,38 +2,36 @@ import numpy as np import pandas as pd -import xarray as xr from huggingface_hub import HfFileSystem fs = HfFileSystem() def get_pv_metadata(testset: pd.DataFrame): - # download from hugginface or load from cache + # download from huggingface or load from cache cache_dir = "data/pv" metadata_file = f"{cache_dir}/metadata.csv" + if not os.path.exists(metadata_file): os.makedirs(cache_dir, exist_ok=True) fs.get("datasets/openclimatefix/uk_pv/metadata.csv", metadata_file) # Load in the dataset metadata_df = pd.read_csv(metadata_file) + metadata_df = metadata_df.rename(columns={"ss_id": "pv_id"}) # join metadata with testset - metadata_df = metadata_df.rename(columns={"ss_id": "pv_id"}) combined_data = testset.merge(metadata_df, on="pv_id", how="left") - # only keep the columns we need + # Select and rename columns combined_data = combined_data[ - ["pv_id", "timestamp", "latitude_rounded", "longitude_rounded", "kwp"] + ["pv_id", "timestamp", "latitude_rounded", "longitude_rounded", "kWp"] ] - - # rename latitude_rounded to latitude and longitude_rounded to longitude combined_data = combined_data.rename( columns={ "latitude_rounded": "latitude", "longitude_rounded": "longitude", - "kwp": "capacity", + "kWp": "capacity", } ) @@ -44,50 +42,91 @@ def get_pv_metadata(testset: pd.DataFrame): def get_pv_truth(testset: pd.DataFrame): + """ + Load PV ground truth data from Hugging Face dataset. + Dataset uses parquet format: 5_minutely/year=YYYY/month=MM/data.parquet + """ print("Loading PV data") - # download from hugginface or load from cache - cache_dir = "data/pv" - metadata_file = f"{cache_dir}/pv.netcdf" - if not os.path.exists(metadata_file): - print("Loading from HF)") - os.makedirs(cache_dir, exist_ok=True) - fs.get("datasets/openclimatefix/uk_pv/pv.netcdf", metadata_file) + cache_dir = "data/pv/parquet_cache" + os.makedirs(cache_dir, exist_ok=True) + + # Get unique year-month combinations from testset + testset["timestamp_dt"] = pd.to_datetime(testset["timestamp"]) + testset["year"] = testset["timestamp_dt"].dt.year + testset["month"] = testset["timestamp_dt"].dt.month + year_months = testset[["year", "month"]].drop_duplicates() + + # downloads and load parquet files for each year-month + all_pv_data = [] + for _, row in year_months.iterrows(): + year = int(row["year"]) + month = int(row["month"]) + cache_file = f"{cache_dir}/pv_{year}_{month:02d}.parquet" + + if not os.path.exists(cache_file): + print(f"Downloading {year}-{month:02d}") + hf_path = ( + f"datasets/openclimatefix/uk_pv/5_minutely" + f"/year={year}/month={month:02d}/data.parquet" + ) + fs.get(hf_path, cache_file) - # Load in the dataset - pv_ds = xr.open_dataset(metadata_file, engine="h5netcdf") + df = pd.read_parquet(cache_file) + all_pv_data.append(df) + print(f"Loaded {len(df)} records from {year}-{month:02d}") + # Combine and prepare data + pv_data = pd.concat(all_pv_data, ignore_index=True) + pv_data = pv_data.rename( + columns={ + "ss_id": "pv_id", + "datetime_GMT": "timestamp", + "generation_Wh": "generation_wh", + } + ) + pv_data["timestamp"] = pd.to_datetime(pv_data["timestamp"]) + pv_data["value"] = pv_data["generation_wh"] / 1000 # Convert Wh to kW + + # Generate forecast horizons for each testset entry combined_data = [] for index, row in testset.iterrows(): - print(f"Processing {index} of {len(testset)}") - pv_id = str(row["pv_id"]) + print(f"Processing {index + 1} of {len(testset)}") + pv_id = row["pv_id"] base_datetime = pd.to_datetime(row["timestamp"]) - # Calculate future timestamps up to the max horizon - for i in range(0, 49): # 48 hours in steps of 1 hour - future_datetime = base_datetime + pd.DateOffset(hours=i) - horizon = i # horizon in hours - - try: - # Attempt to select data for the future datetime - selected_data = pv_ds[pv_id].sel(datetime=future_datetime) - value = selected_data.values.item() - value = value / 1000 # to convert from w to kw - except KeyError: - # If data is not found for the future datetime, set value as NaN + # Match timezone with PV data + if base_datetime.tz is None and pv_data["timestamp"].dt.tz is not None: + base_datetime = base_datetime.tz_localize("UTC") + + # Generate 48-hour forecast horizon (0-48 hours) + for i in range(49): + future_datetime = base_datetime + pd.Timedelta(hours=i) + time_window = pd.Timedelta(minutes=5) + + # finds closest matching timestamp within 5-minute window + mask = ( + (pv_data["pv_id"] == pv_id) + & (pv_data["timestamp"] >= future_datetime - time_window) + & (pv_data["timestamp"] <= future_datetime + time_window) + ) + matching_data = pv_data[mask] + + if len(matching_data) > 0: + # find closest match by time difference + matching_data = matching_data.copy() + matching_data["time_diff"] = abs(matching_data["timestamp"] - future_datetime) + value = matching_data.loc[matching_data["time_diff"].idxmin(), "value"] + else: value = np.nan - # Add the data to the DataFrame combined_data.append( - pd.DataFrame( - { - "pv_id": pv_id, - "timestamp": future_datetime, - "value": value, - "horizon_hour": horizon, - }, - index=[i], - ) + { + "pv_id": pv_id, + "timestamp": future_datetime, + "value": value, + "horizon_hour": i, + } ) - combined_data = pd.concat(combined_data) - return combined_data + + return pd.DataFrame(combined_data) diff --git a/quartz_solar_forecast/evaluation.py b/quartz_solar_forecast/evaluation.py index ded0843af..f04439a1d 100644 --- a/quartz_solar_forecast/evaluation.py +++ b/quartz_solar_forecast/evaluation.py @@ -31,7 +31,7 @@ ) -def run_eval(testset_path: str = "dataset/testset.csv"): +def run_eval(testset_path: str = os.path.join(os.path.dirname(__file__), "dataset", "testset.csv")): # load testset from csv testset = pd.read_csv(testset_path) @@ -63,4 +63,4 @@ def run_eval(testset_path: str = "dataset/testset.csv"): # TODO -# run_eval() +run_eval()