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")