forecast/forecast-algorithm/短期中期部署/部署/OuLiPowerV1.0/utils/weather_utils.py

import math

import numpy as np
import pandas as pd
from typing import Tuple

def calculate_wind_speed_direction(u: float, v: float) -> Tuple[float, float]:
    """
    Calculate wind speed and direction from u and v wind components.
    
    Args:
        u (float): East-west wind component (positive for eastward wind)
        v (float): North-south wind component (positive for northward wind)
        
    Returns:
        Tuple[float, float]: Wind speed (m/s) and direction (degrees from north, clockwise)
        
    Note:
        - Wind direction is in meteorological convention (degrees from north, clockwise)
        - Wind direction indicates where the wind is coming FROM (not going to)
        - 0 degrees is from the north
        - 90 degrees is from the east
        - 180 degrees is from the south
        - 270 degrees is from the west
    """
    # Calculate wind speed
    wind_speed = np.sqrt(u**2 + v**2)

    # Calculate wind direction using meteorological convention
    # arctan2(x, y) gives the angle in radians between the positive x-axis and the point (x, y)
    # For wind direction (where wind is coming FROM), we need to use negative components
    # and adjust the angle to get the meteorological convention
    wind_direction_rad = np.arctan2(-u, -v)
    wind_direction_deg = np.degrees(wind_direction_rad)
    wind_direction_deg = (wind_direction_deg + 360) % 360
    
    return wind_speed, wind_direction_deg

def calculate_wind_components(speed: float, direction: float) -> Tuple[float, float]:
    """
    Calculate u and v wind components from wind speed and direction.
    
    Args:
        speed (float): Wind speed in m/s
        direction (float): Wind direction in degrees from north, clockwise (where wind is coming FROM)
        
    Returns:
        Tuple[float, float]: u (east-west) and v (north-south) wind components
        
    Note:
        - Positive u means eastward wind
        - Positive v means northward wind
        - This is the inverse operation of calculate_wind_speed_direction
    """
    # Convert direction from degrees to radians
    # Since direction is where wind is coming FROM, we need to add 180 degrees to get where it's going TO
    direction_rad = np.radians((direction + 180) % 360)
    
    # Calculate components
    u = -speed * np.sin(direction_rad)
    v = -speed * np.cos(direction_rad)
    
    return u, v 


def interpolate_weather_attributes(df: pd.DataFrame, time_col: str, time_interval: str='15min', method: str='cubic') -> pd.DataFrame:
    """
    Interpolate weather attributes to a uniform time interval using specified method.
    
    This function handles cases where input data has varying frequencies 
    (e.g., first 72h hourly, next 72h 3-hourly, last 72h 6-hourly) and 
    creates a uniform time grid with interpolation.
    
    Args:
        df (pd.DataFrame): Input DataFrame with weather attributes
        time_col (str): Name of the time column
        time_interval (str): Target time interval for interpolation (e.g., '15min', '30min', '1h')
        method (str): Interpolation method to use ('linear' or 'cubic')
        
    Returns:
        pd.DataFrame: DataFrame with interpolated weather attributes at uniform intervals
        
    Example:
        df_interp = interpolate_weather_attributes(df, 'time', '15min', method='cubic')
    """
    if df.empty:
        return df.copy()
    
    if method not in ['linear', 'cubic']:
        raise ValueError("Method must be either 'linear' or 'cubic'")
    
    # Make a copy to avoid modifying the original DataFrame
    df_work = df.copy()
    
    # Ensure time column is datetime
    if not pd.api.types.is_datetime64_any_dtype(df_work[time_col]):
        df_work[time_col] = pd.to_datetime(df_work[time_col])
    
    # Sort by time to ensure proper interpolation
    df_work = df_work.sort_values(time_col).reset_index(drop=True)
    
    # Remove any duplicate timestamps, keeping the first occurrence
    df_work = df_work.drop_duplicates(subset=[time_col], keep='first')
    
    # Get the time range
    start_time = df_work[time_col].min()
    end_time = df_work[time_col].max()
    
    # Create uniform time grid
    time_grid = pd.date_range(start=start_time, end=end_time, freq=time_interval)
    
    # Create new DataFrame with uniform time grid
    df_uniform = pd.DataFrame({time_col: time_grid})
    
    # Set time as index for both DataFrames to facilitate interpolation
    df_work_indexed = df_work.set_index(time_col)
    df_uniform_indexed = df_uniform.set_index(time_col)
    
    # Get all columns except the time column for interpolation
    value_columns = [col for col in df_work.columns if col != time_col]
    
    # Perform interpolation for each weather attribute
    for col in value_columns:
        if pd.api.types.is_numeric_dtype(df_work_indexed[col]):
            # For numeric columns, use specified interpolation method
            # First, reindex to include all target time points
            series_reindexed = df_work_indexed[col].reindex(
                df_work_indexed.index.union(df_uniform_indexed.index)
            ).sort_index()
            
            # Interpolate using specified method
            series_interpolated = series_reindexed.interpolate(method=method)
            
            # Extract values at uniform time grid points
            df_uniform[col] = series_interpolated.reindex(df_uniform_indexed.index).values
            
        else:
            # For non-numeric columns, use forward fill
            series_reindexed = df_work_indexed[col].reindex(
                df_work_indexed.index.union(df_uniform_indexed.index)
            ).sort_index()
            
            series_filled = series_reindexed.ffill()
            df_uniform[col] = series_filled.reindex(df_uniform_indexed.index).values
    
    # Handle any remaining NaN values at the beginning or end
    # Forward fill first, then backward fill
    df_uniform = df_uniform.ffill().bfill()
    
    return df_uniform

def nearest_position(stn_lat: float, stn_lon: float, lat2d: np.ndarray, lon2d: np.ndarray) -> Tuple[int, int]:
    """
    Find the nearest grid point to a given station location using great circle distance.
    
    Args:
        stn_lat (float): Station latitude in degrees
        stn_lon (float): Station longitude in degrees
        lat2d (np.ndarray): 2D array of grid latitudes
        lon2d (np.ndarray): 2D array of grid longitudes
        
    Returns:
        Tuple[int, int]: (row, col) indices of the nearest grid point
        
    Note:
        Uses great circle distance calculation which is more accurate for geographic coordinates
        than simple Euclidean distance.
    """
    # Convert to radians for distance calculation
    stn_lat_rad = np.radians(stn_lat)
    stn_lon_rad = np.radians(stn_lon)
    lat2d_rad = np.radians(lat2d)
    lon2d_rad = np.radians(lon2d)
    
    # Calculate great circle distance using haversine formula
    dlat = lat2d_rad - stn_lat_rad
    dlon = lon2d_rad - stn_lon_rad
    a = np.sin(dlat/2)**2 + np.cos(stn_lat_rad) * np.cos(lat2d_rad) * np.sin(dlon/2)**2
    c = 2 * np.arcsin(np.sqrt(a))
    distance = 6371 * c  # Earth's radius in km
    
    # Find indices of minimum distance
    min_idx = np.unravel_index(np.argmin(distance), distance.shape)
    
    return min_idx


def calculate_sunrise_sunset(date, latitude, longitude):
    """
    Calculate sunrise and sunset times for a given date and location.
    Based on the equations from NOAA:
    https://www.esrl.noaa.gov/gmd/grad/solcalc/

    Args:
        date: Date to calculate sunrise/sunset for (in Beijing time)
        latitude: Latitude in degrees
        longitude: Longitude in degrees

    Returns:
        Tuple of (sunrise_hour, sunset_hour) in Beijing time (decimal hours)
    """
    # Convert latitude and longitude to radians
    lat_rad = math.radians(latitude)

    # Calculate day of year
    day_of_year = date.timetuple().tm_yday

    # Calculate solar declination angle (in radians)
    declination = 0.4093 * math.sin(2 * math.pi / 365 * (day_of_year - 81))

    # Calculate equation of time (in minutes)
    B = 2 * math.pi / 365 * (day_of_year - 81)
    time_eq = 9.87 * math.sin(2 * B) - 7.53 * math.cos(B) - 1.5 * math.sin(B)

    # Calculate solar hour angle at sunrise/sunset (in radians)
    # Using -0.83 degrees as the solar elevation angle at sunrise/sunset
    solar_elevation = math.radians(-0.83)
    ha = math.acos(
        (math.sin(solar_elevation) - math.sin(lat_rad) * math.sin(declination)) /
        (math.cos(lat_rad) * math.cos(declination))
    )

    # Convert hour angle to hours
    ha_hours = math.degrees(ha) / 15.0

    # Calculate sunrise and sunset (in local solar time, hours)
    sunrise_local = 12.0 - ha_hours - time_eq / 60.0
    sunset_local = 12.0 + ha_hours - time_eq / 60.0

    # Convert local solar time to Beijing time (UTC+8, 120°E)
    # Beijing time = local solar time + (120 - longitude)/15
    longitude_correction = (120.0 - longitude) / 15.0
    sunrise_beijing = (sunrise_local + longitude_correction) % 24
    sunset_beijing = (sunset_local + longitude_correction) % 24

    return sunrise_beijing, sunset_beijing


def solar_elevation_angle(date_time, latitude, longitude, timezone_longitude=120):
    """
    Calculate solar elevation angle in degrees.

    Args:
        date_time: datetime object
        latitude: latitude in degrees (-90 to 90)
        longitude: longitude in degrees (-180 to 180)
        timezone_longitude: timezone longitude in degrees (-180 to 180)

    Returns:
        Solar elevation angle in degrees
    """
    # Input validation
    if not -90 <= latitude <= 90:
        raise ValueError("Latitude must be between -90 and 90 degrees")
    if not -180 <= longitude <= 180:
        raise ValueError("Longitude must be between -180 and 180 degrees")
    if not -180 <= timezone_longitude <= 180:
        raise ValueError("Timezone longitude must be between -180 and 180 degrees")

    # Time information
    n = date_time.timetuple().tm_yday
    decimal_hours = date_time.hour + date_time.minute / 60 + date_time.second / 3600

    # Calculate solar declination and hour angle
    delta = solar_declination(n)
    H = hour_angle(decimal_hours, longitude, timezone_longitude, n)

    # Calculate elevation angle
    phi = latitude
    sin_h = math.sin(math.radians(phi)) * math.sin(math.radians(delta)) + \
            math.cos(math.radians(phi)) * math.cos(math.radians(delta)) * math.cos(math.radians(H))
    h = math.degrees(math.asin(sin_h))

    # Apply atmospheric refraction correction
    h += atmospheric_refraction(h)

    return h


def solar_declination(n):
    """Calculate solar declination angle in degrees"""
    return 23.45 * math.sin(math.radians(360 / 365 * (n - 81)))


def hour_angle(local_time, longitude, timezone_longitude, n):
    """Calculate hour angle in degrees"""
    eot = equation_of_time(n)
    lst = local_time + (timezone_longitude - longitude) * 4 / 60 + eot / 60
    return 15 * (lst - 12)


def equation_of_time(n):
    """Calculate the equation of time in minutes"""
    B = 360 / 365 * (n - 81)
    return 9.87 * math.sin(2 * math.radians(B)) - 7.53 * math.cos(math.radians(B)) - 1.5 * math.sin(math.radians(B))


def atmospheric_refraction(elevation):
    """
    Calculate atmospheric refraction correction in degrees using Bennett's formula.

    Atmospheric refraction causes celestial objects to appear higher in the sky
    than their true geometric position. This effect is strongest at the horizon
    and decreases with increasing elevation angle.

    Args:
        elevation: Apparent elevation angle in degrees

    Returns:
        Atmospheric refraction correction in degrees (always positive)

    Reference:
        Bennett, G.G. (1982). "The Calculation of Astronomical Refraction in Marine Navigation"
        Journal of Navigation, 35(2), 255-259
    """
    # No refraction for objects below the horizon (accounting for refraction itself)
    if elevation < -0.575:
        return 0

    # For very high elevations, refraction is negligible
    if elevation > 85:
        return 0

    # Bennett's formula for atmospheric refraction
    # R = 1.02 / tan(h + 10.3/(h + 5.11)) arcminutes
    # Convert to degrees by dividing by 60

    # Handle low elevations with polynomial approximation for better accuracy
    if elevation < 15:
        # Sæmundsson's formula for low elevations (more accurate near horizon)
        # R in arcseconds = 1735 + h*(-518.2 + h*(103.4 + h*(-12.79 + h*0.711)))
        # where h is in degrees
        h = elevation
        refraction_arcsec = 1735 + h * (-518.2 + h * (103.4 + h * (-12.79 + h * 0.711)))
        refraction_deg = refraction_arcsec / 3600  # Convert arcseconds to degrees
    else:
        # Bennett's formula for higher elevations
        h = elevation
        refraction_arcmin = 1.02 / math.tan(math.radians(h + 10.3 / (h + 5.11)))
        refraction_deg = refraction_arcmin / 60  # Convert arcminutes to degrees

    return refraction_deg


def solar_azimuth_angle(date_time, latitude, longitude, timezone_longitude=120):
    """
    Calculate solar azimuth angle in degrees.

    The azimuth angle is measured clockwise from true North:
    - 0° = North
    - 90° = East
    - 180° = South
    - 270° = West

    Args:
        date_time: datetime object
        latitude: latitude in degrees (-90 to 90)
        longitude: longitude in degrees (-180 to 180)
        timezone_longitude: timezone longitude in degrees (-180 to 180)

    Returns:
        Solar azimuth angle in degrees (0-360)

    Reference:
        Uses standard solar position equations from NOAA and astronomical texts
    """
    # Input validation
    if not -90 <= latitude <= 90:
        raise ValueError("Latitude must be between -90 and 90 degrees")
    if not -180 <= longitude <= 180:
        raise ValueError("Longitude must be between -180 and 180 degrees")
    if not -180 <= timezone_longitude <= 180:
        raise ValueError("Timezone longitude must be between -180 and 180 degrees")

    # Time information
    n = date_time.timetuple().tm_yday
    decimal_hours = date_time.hour + date_time.minute / 60 + date_time.second / 3600

    # Calculate solar declination and hour angle
    delta = solar_declination(n)  # Solar declination in degrees
    H = hour_angle(decimal_hours, longitude, timezone_longitude, n)  # Hour angle in degrees

    # Convert to radians for calculation
    phi_rad = math.radians(latitude)
    delta_rad = math.radians(delta)
    H_rad = math.radians(H)

    # Calculate azimuth using atan2 for correct quadrant handling
    # Formula: A = atan2(sin(H), cos(H) × sin(φ) - tan(δ) × cos(φ))
    # This gives azimuth measured from South, positive towards West

    # Calculate elevation angle for reference
    sin_h = math.sin(phi_rad) * math.sin(delta_rad) + \
            math.cos(phi_rad) * math.cos(delta_rad) * math.cos(H_rad)

    # Calculate azimuth components
    numerator = math.sin(H_rad)
    denominator = math.cos(H_rad) * math.sin(phi_rad) - math.tan(delta_rad) * math.cos(phi_rad)

    # Calculate azimuth from South
    azimuth_rad = math.atan2(numerator, denominator)
    azimuth_deg = math.degrees(azimuth_rad)

    # Convert from "from South" to "from North" convention
    # Azimuth from South: 0° = South, positive towards West
    # We want: 0° = North, positive clockwise (East)
    azimuth_from_north = (azimuth_deg + 180) % 360

    return azimuth_from_north


def calculate_clear_sky_radiation(date_time, latitude, longitude, timezone_longitude=120):
    """
    Calculate theoretical clear-sky solar radiation using simplified Bird Clear Sky Model.

    This represents the maximum possible radiation for given sun position under
    perfectly clear atmospheric conditions.

    Args:
        date_time: datetime object
        latitude: latitude in degrees (-90 to 90)
        longitude: longitude in degrees (-180 to 180)
        timezone_longitude: timezone longitude in degrees (-180 to 180)

    Returns:
        Clear-sky radiation in W/m² (0 if sun is below horizon)

    Reference:
        Simplified Bird Clear Sky Model
        R.E. Bird and R.L. Hulstrom (1981)
    """
    # Calculate solar elevation angle (already includes atmospheric refraction)
    elevation = solar_elevation_angle(date_time, latitude, longitude, timezone_longitude)

    # If sun is below horizon, no radiation
    if elevation <= 0:
        return 0.0

    # Solar constant (W/m²)
    solar_constant = 1367.0

    # Day of year for seasonal correction
    n = date_time.timetuple().tm_yday

    # Earth-Sun distance correction factor
    # Accounts for elliptical orbit
    B = 2 * math.pi * (n - 1) / 365
    distance_factor = 1.000110 + 0.034221 * math.cos(B) + 0.001280 * math.sin(B) + \
                     0.000719 * math.cos(2*B) + 0.000077 * math.sin(2*B)

    # Extraterrestrial radiation on horizontal surface
    elevation_rad = math.radians(elevation)
    I0 = solar_constant * distance_factor * math.sin(elevation_rad)

    # Simplified atmospheric transmittance model
    # Air mass (Kasten and Young formula)
    air_mass = 1.0 / (math.sin(elevation_rad) + 0.50572 * (elevation + 6.07995)**(-1.6364))

    # Simplified atmospheric transmittance (assumes standard atmosphere)
    # Typical clear-sky conditions: ~0.7-0.8 transmittance
    # Components: Rayleigh scattering, aerosol extinction, water vapor, ozone
    tau_rayleigh = math.exp(-0.0903 * air_mass**0.84)
    tau_aerosol = math.exp(-0.0127 * air_mass**0.26)  # Clean atmosphere assumption
    tau_water = math.exp(-0.0365 * air_mass**0.5)     # Standard water vapor

    # Total transmittance
    transmittance = tau_rayleigh * tau_aerosol * tau_water

    # Clear-sky direct + diffuse radiation
    # Direct beam
    direct_beam = I0 * transmittance

    # Diffuse radiation (simplified - roughly 10-20% of direct for clear sky)
    diffuse_factor = 0.15
    diffuse = I0 * diffuse_factor

    # Total clear-sky radiation
    clear_sky_ghi = direct_beam + diffuse

    # Ensure non-negative
    return max(0.0, clear_sky_ghi)


if __name__ == "__main__":
    # Test wind speed/direction calculations
    # u = 1.3415451049804688  # m/s eastward
    # v = -2.3990478515625 # m/s northward
    # speed, direction = calculate_wind_speed_direction(u, v)
    # print(f"Wind speed: {speed:.2f} m/s")
    # print(f"Wind direction: {direction:.2f}°")
    
    # Test nearest_position function
    # print("\n=== Testing Nearest Position Function ===")
    # # Create a sample grid
    # lats = np.linspace(45.0, 46.0, 10)  # 10 points from 45°N to 46°N
    # lons = np.linspace(84.0, 85.0, 10)  # 10 points from 84°E to 85°E
    # lat2d, lon2d = np.meshgrid(lats, lons)
    
    # # Test point (close to the target location in the notebook)
    # test_lat = 45.94375072
    # test_lon = 84.22134241
    
    # # Find nearest point
    # row, col = nearest_position(test_lat, test_lon, lat2d, lon2d)
    # print(f"Test point: ({test_lat}, {test_lon})")
    # print(f"Nearest grid point: ({lat2d[row, col]:.6f}, {lon2d[row, col]:.6f})")
    # print(f"Grid indices: row={row}, col={col}")
    
    # Test interpolation function with varying time frequencies
    print("\n=== Testing Weather Data Interpolation ===")
    
    # Create sample data with varying frequencies
    # First 6 hours: hourly data
    time1 = pd.date_range('2023-01-01 00:00', periods=7, freq='1h')
    temp1 = np.linspace(15, 18, 7)  # Temperature rising from 15 to 18°C
    ws1 = np.linspace(5, 8, 7)      # Wind speed increasing from 5 to 8 m/s
    
    # Next 6 hours: 3-hourly data  
    time2 = pd.date_range('2023-01-01 06:00', periods=3, freq='3h')
    temp2 = np.linspace(18, 22, 3)  # Temperature rising from 18 to 22°C
    ws2 = np.linspace(8, 12, 3)     # Wind speed increasing from 8 to 12 m/s
    
    # Last 6 hours: 6-hourly data
    time3 = pd.date_range('2023-01-01 12:00', periods=2, freq='6h')
    temp3 = np.linspace(22, 20, 2)  # Temperature dropping from 22 to 20°C
    ws3 = np.linspace(12, 10, 2)    # Wind speed decreasing from 12 to 10 m/s
    
    # Combine all data
    times = np.concatenate([time1, time2[1:], time3[1:]])  # Remove overlapping points
    temperatures = np.concatenate([temp1, temp2[1:], temp3[1:]])
    wind_speeds = np.concatenate([ws1, ws2[1:], ws3[1:]])
    
    # Create DataFrame
    df_varying = pd.DataFrame({
        'time': times,
        'temperature': temperatures, 
        'wind_speed': wind_speeds,
    })
    
    print("Original data with varying frequencies:")
    print(df_varying)
    print(f"\nOriginal data shape: {df_varying.shape}")
    
    # Interpolate to 30-minute intervals
    df_interpolated = interpolate_weather_attributes(df_varying, 'time', '15min', method='linear')
    
    print(f"\nInterpolated data (30-min intervals):")
    print(df_interpolated.head(100))
    print(f"Interpolated data shape: {df_interpolated.shape}")
    
    print(f"\nTime range: {df_interpolated['time'].min()} to {df_interpolated['time'].max()}")
    print(f"Number of time points: {len(df_interpolated)}")
    
    # # u, v = calculate_wind_components(speed, direction)
    # # print(f"u component: {u:.2f} m/s")
    # # print(f"v component: {v:.2f} m/s")