pvlib中累积GHI的时间约定问题咨询

Great question—this is such a common pitfall when working with hourly aggregated meteorological data and pvlib's strict focus on instantaneous (or mid-hour representative) time steps. Let's break down exactly how to fix the time offset issue and use cumulative GHI properly.

The Core Problem #

pvlib's functions (like pvlib.irradiance.erbs) expect irradiance values (W/m²) tied to a specific moment in time (or a representative mid-point for hourly data). GFS's hourly GHI, however, is a cumulative energy value (J/m²) for the previous hour (e.g., the 12:00 timestamp represents 11:00–12:00 total radiation). Using the instantaneous zenith at 12:00 with this cumulative GHI misaligns the solar position with the actual period the radiation was measured over—hence the time offset in your results.

Step-by-Step Solution #

Here's how to align everything correctly:

1. Convert Cumulative GHI to Average Hourly Irradiance
   First, turn the cumulative energy (J/m²) into average irradiance (W/m²) by dividing by 3600 seconds (the length of the hour):

   df['ghi_avg'] = df['ghi'] / 3600
   

   This gives you the average power per square meter over the hour, which is the unit pvlib expects.

2. Shift Time Index to the Hour Midpoint
   Move your timestamp from the top of the hour (e.g., 12:00) to the middle (12:30). This ensures the solar position you calculate corresponds to the average conditions of the hour, matching the window your cumulative GHI covers:

   df_mid = df.copy()
   df_mid.index = df_mid.index + pd.Timedelta(minutes=30)
   

3. Recalculate Solar Position for Mid-Hour Timestamps
   Use the shifted index to get solar position values that represent the hour's average conditions. Always use apparent zenith for irradiance calculations—it accounts for atmospheric refraction:

   solpos = pvlib.solarposition.get_solarposition(
       time=df_mid.index,
       latitude=your_lat,
       longitude=your_lon
   )
   

4. Run pvlib Functions with Aligned Data
   Now feed the mid-hour irradiance and solar position into erbs (or any other pvlib function):

   df_dni_dhi_kt = pvlib.irradiance.erbs(
       ghi=df_mid['ghi_avg'],
       zenith=solpos['apparent_zenith'],
       datetime_or_doy=df_mid.index
   )
   

5. Optional: Align Back to Hourly Timestamps
   If you need to compare your results to hourly observations tied to top-of-hour timestamps, shift the index back:

   df_dni_dhi_kt.index = df_dni_dhi_kt.index - pd.Timedelta(minutes=30)
   

Key Notes #

• Verify GFS's Time Window: Double-check that GFS's hourly GHI is indeed the previous hour (not future). Some models label hourly data as the end of the window, others as the start—this affects how you shift timestamps.
• Mid-Hour vs. Average Solar Position: For most use cases, mid-hour solar position is accurate enough. If you need extreme precision, you could calculate solar position at the start and end of the hour and take the average—but mid-hour is simpler and works for 99% of applications.
• Power Calculation Alignment: When you get to AC/DC power calculations, ensure your load/observation data uses the same time convention (mid-hour or shifted back to top-of-hour) to avoid offsets.

Full Example Code #

import pandas as pd
import pvlib

# Sample GFS data: top-of-hour timestamps, cumulative GHI (J/m²)
gfs_df = pd.DataFrame(
    {'ghi': [216000, 432000, 648000]},  # 60W/m², 120W/m², 180W/m² hourly average
    index=pd.date_range('2024-06-21 09:00', periods=3, freq='H')
)
latitude = 35.7796
longitude = -78.6382

# 1. Convert cumulative to average irradiance
gfs_df['ghi_avg'] = gfs_df['ghi'] / 3600

# 2. Shift to mid-hour timestamps
gfs_mid = gfs_df.copy()
gfs_mid.index = gfs_mid.index + pd.Timedelta(minutes=30)

# 3. Get aligned solar position
solpos = pvlib.solarposition.get_solarposition(
    time=gfs_mid.index,
    latitude=latitude,
    longitude=longitude
)

# 4. Run erbs with aligned data
irrad_components = pvlib.irradiance.erbs(
    ghi=gfs_mid['ghi_avg'],
    zenith=solpos['apparent_zenith'],
    datetime_or_doy=gfs_mid.index
)

# 5. Shift back to top-of-hour if needed
irrad_components.index = irrad_components.index - pd.Timedelta(minutes=30)

print(irrad_components)

内容的提问来源于stack exchange，提问作者JulienV