Spatial Interpolation

def natural_neighbor_to_grid(xp, yp, variable, interp_type='linear', hres=50000, 
                           search_radius=None, rbf_func='linear', rbf_smooth=0,
                           minimum_neighbors=3, gamma=0.25, kappa_star=5.052, 
                           boundary_coords=None, grid_projection=None):
    """
    Generate a natural neighbor interpolation of the given point data to a grid.
    
    Natural neighbor interpolation is a method of spatial interpolation for
    multivariate data that preserves local properties and is well-suited for
    irregular station networks common in meteorology.
    
    Parameters:
    - xp: x-coordinates of observation points
    - yp: y-coordinates of observation points  
    - variable: data values at observation points
    - interp_type: interpolation method ('linear', 'nearest', 'cubic')
    - hres: horizontal resolution of output grid
    - search_radius: maximum search radius for points
    - rbf_func: radial basis function type
    - rbf_smooth: smoothing parameter for RBF
    - minimum_neighbors: minimum neighbors required for interpolation
    - gamma: shape parameter for natural neighbor weighting
    - kappa_star: scaling parameter for influence radius
    - boundary_coords: boundary coordinates for interpolation domain
    - grid_projection: coordinate reference system for output grid
    
    Returns:
    Tuple of (grid_x, grid_y, grid_data) arrays with interpolated values
    """

def inverse_distance_to_grid(xp, yp, variable, hres, search_radius=None, 
                           rbf_func='linear', rbf_smooth=0, minimum_neighbors=3,
                           gamma=0.25, kappa_star=5.052, boundary_coords=None):
    """
    Interpolate point data to grid using inverse distance weighting.
    
    Inverse distance weighting is a simple but effective interpolation method
    that weights observations inversely proportional to their distance from
    the interpolation point.
    
    Parameters:
    - xp: x-coordinates of observation points
    - yp: y-coordinates of observation points
    - variable: data values at observation points
    - hres: horizontal resolution of output grid
    - search_radius: maximum search radius for points
    - rbf_func: radial basis function type
    - rbf_smooth: smoothing parameter
    - minimum_neighbors: minimum neighbors for interpolation
    - gamma: weighting exponent (typically 1-3)
    - kappa_star: influence radius scaling
    - boundary_coords: interpolation domain boundary
    
    Returns:
    Tuple of (grid_x, grid_y, grid_data) arrays
    """

def barnes_to_grid(xp, yp, variable, hres, search_radius, kappa_star=5.052,
                  gamma=0.25, minimum_neighbors=3, boundary_coords=None):
    """
    Interpolate point data using Barnes objective analysis.
    
    Barnes analysis is a successive correction method widely used in meteorology
    for creating gridded analyses from irregular station data. It applies
    Gaussian weighting with multiple passes to reduce analysis error.
    
    Parameters:
    - xp: x-coordinates of observation points
    - yp: y-coordinates of observation points
    - variable: data values at observation points
    - hres: horizontal resolution of output grid
    - search_radius: search radius for station selection
    - kappa_star: Gaussian response parameter (typically 3-10)
    - gamma: convergence parameter for successive corrections
    - minimum_neighbors: minimum stations required for analysis
    - boundary_coords: analysis domain boundary coordinates
    
    Returns:
    Tuple of (grid_x, grid_y, grid_data) with Barnes analysis
    """

def cressman_to_grid(xp, yp, variable, hres, search_radius, minimum_neighbors=3,
                    boundary_coords=None):
    """
    Interpolate point data using Cressman objective analysis.
    
    Cressman analysis uses successive corrections with distance-weighted
    influence functions. It's computationally efficient and widely used
    in operational meteorology for real-time analysis.
    
    Parameters:
    - xp: x-coordinates of observation points
    - yp: y-coordinates of observation points
    - variable: data values at observation points
    - hres: horizontal resolution of output grid
    - search_radius: influence radius for each correction pass
    - minimum_neighbors: minimum observations for grid point analysis
    - boundary_coords: analysis domain boundaries
    
    Returns:
    Tuple of (grid_x, grid_y, grid_data) with Cressman analysis
    """

def geodetic_to_grid(xp, yp, variable, boundary_coords, hres, interp_type='linear',
                    search_radius=None, rbf_func='linear', rbf_smooth=0,
                    minimum_neighbors=3, gamma=0.25, grid_projection=None):
    """
    Interpolate data from geodetic coordinates (lat/lon) to a projected grid.
    
    Handles coordinate transformation from geographic to projected coordinate
    systems while performing spatial interpolation, essential for meteorological
    analysis spanning large geographic areas.
    
    Parameters:
    - xp: longitude coordinates of observation points
    - yp: latitude coordinates of observation points
    - variable: data values at observation points
    - boundary_coords: geographic bounds for interpolation domain
    - hres: horizontal resolution in projected coordinates
    - interp_type: interpolation method
    - search_radius: search radius in projected coordinates
    - rbf_func: radial basis function type
    - rbf_smooth: smoothing parameter
    - minimum_neighbors: minimum neighbors for interpolation
    - gamma: method-specific parameter
    - grid_projection: target projection for output grid
    
    Returns:
    Tuple of (grid_x, grid_y, grid_data) in projected coordinates
    """

def remove_nan_observations(x, y, *variables):
    """
    Remove observations with NaN values from coordinate and data arrays.
    
    Quality control step essential for reliable interpolation results,
    removing invalid or missing observations that would contaminate
    the analysis.
    
    Parameters:
    - x: x-coordinate array
    - y: y-coordinate array  
    - variables: one or more data variable arrays
    
    Returns:
    Tuple of cleaned (x, y, *variables) with NaN observations removed
    """

def remove_repeat_coordinates(x, y, z, radius=None):
    """
    Remove observations with duplicate or nearly duplicate coordinates.
    
    Eliminates station observations that are too close together, which
    can cause numerical issues in interpolation algorithms and bias
    the analysis toward over-sampled regions.
    
    Parameters:
    - x: x-coordinate array
    - y: y-coordinate array
    - z: data value array
    - radius: minimum separation distance between stations
    
    Returns:
    Tuple of (x, y, z) with duplicate coordinates removed
    """

def interpolate_1d(x, xp, *variables, axis=0, fill_value=np.nan):
    """
    Interpolate 1D arrays to new coordinate values.
    
    Essential for vertical interpolation of atmospheric soundings to
    standard pressure levels, height coordinates, or potential temperature
    surfaces commonly used in meteorological analysis.
    
    Parameters:
    - x: new coordinate values for interpolation
    - xp: original coordinate values (must be monotonic)
    - variables: one or more data arrays to interpolate
    - axis: axis along which to interpolate
    - fill_value: value for extrapolation outside data range
    
    Returns:
    Interpolated arrays at new coordinate values
    """

def interpolate_nans(array, method='linear', axis=0, limit=None):
    """
    Fill NaN values in array using interpolation.
    
    Useful for filling missing values in time series or vertical profiles
    while preserving the overall structure and variability of the data.
    
    Parameters:
    - array: input array with NaN values to fill
    - method: interpolation method ('linear', 'nearest', 'cubic')
    - axis: axis along which to interpolate
    - limit: maximum number of consecutive NaNs to fill
    
    Returns:
    Array with NaN values filled by interpolation
    """

def cross_section(data, start, end, steps=100, interp_type='linear'):
    """
    Extract a cross-section through gridded 3D atmospheric data.
    
    Creates vertical cross-sections along arbitrary horizontal paths,
    essential for analyzing atmospheric structure, frontal zones,
    and three-dimensional meteorological phenomena.
    
    Parameters:
    - data: 3D xarray DataArray with atmospheric data
    - start: starting point coordinates (lat, lon) or (x, y)
    - end: ending point coordinates (lat, lon) or (x, y)
    - steps: number of points along cross-section path
    - interp_type: interpolation method for extracting values
    
    Returns:
    xarray DataArray containing the cross-section with distance and
    vertical coordinates
    """

import metpy.interpolate as mpinterp
from metpy.units import units
import numpy as np
import matplotlib.pyplot as plt

# Sample station data
station_lons = np.array([-100, -95, -90, -85, -80]) * units.degrees_east
station_lats = np.array([35, 40, 45, 50, 45]) * units.degrees_north  
temperatures = np.array([25, 20, 15, 10, 12]) * units.celsius

# Convert to projected coordinates (example: meters)
x_coords = station_lons.m * 111320  # Rough conversion for demo
y_coords = station_lats.m * 111320

# Natural neighbor interpolation to regular grid
grid_x, grid_y, temp_grid = mpinterp.natural_neighbor_to_grid(
    x_coords, y_coords, temperatures.m, 
    hres=50000,  # 50 km resolution
    interp_type='linear'
)

# Plot results
plt.figure(figsize=(10, 8))
plt.contourf(grid_x, grid_y, temp_grid, levels=15, cmap='RdYlBu_r')
plt.colorbar(label='Temperature (°C)')
plt.scatter(x_coords, y_coords, c='black', s=50, label='Stations')
plt.xlabel('X (m)')
plt.ylabel('Y (m)')
plt.title('Temperature Analysis - Natural Neighbor Interpolation')
plt.legend()
plt.show()

# More sophisticated analysis with Barnes method
search_radius = 200000  # 200 km search radius

grid_x, grid_y, temp_barnes = mpinterp.barnes_to_grid(
    x_coords, y_coords, temperatures.m,
    hres=25000,  # 25 km resolution
    search_radius=search_radius,
    kappa_star=5.052,  # Standard Barnes parameter
    gamma=0.25,        # Convergence parameter
    minimum_neighbors=3
)

# Compare with Cressman analysis
grid_x_cress, grid_y_cress, temp_cressman = mpinterp.cressman_to_grid(
    x_coords, y_coords, temperatures.m,
    hres=25000,
    search_radius=search_radius,
    minimum_neighbors=3
)

# Plot comparison
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))

# Barnes analysis
cs1 = ax1.contourf(grid_x, grid_y, temp_barnes, levels=15, cmap='RdYlBu_r')
ax1.scatter(x_coords, y_coords, c='black', s=50)
ax1.set_title('Barnes Analysis')
plt.colorbar(cs1, ax=ax1, label='Temperature (°C)')

# Cressman analysis  
cs2 = ax2.contourf(grid_x_cress, grid_y_cress, temp_cressman, levels=15, cmap='RdYlBu_r')
ax2.scatter(x_coords, y_coords, c='black', s=50)
ax2.set_title('Cressman Analysis')
plt.colorbar(cs2, ax=ax2, label='Temperature (°C)')

plt.tight_layout()
plt.show()

# Example with noisy data requiring quality control
station_x = np.array([100, 150, 200, 200.1, 250, 300]) * 1000  # Some duplicates
station_y = np.array([200, 250, 300, 300.05, 350, 400]) * 1000
station_data = np.array([15.2, 18.7, np.nan, 22.1, 19.8, 16.5])  # Some NaNs

# Remove NaN observations
clean_x, clean_y, clean_data = mpinterp.remove_nan_observations(
    station_x, station_y, station_data
)
print(f"Removed {len(station_data) - len(clean_data)} NaN observations")

# Remove duplicate coordinates
final_x, final_y, final_data = mpinterp.remove_repeat_coordinates(
    clean_x, clean_y, clean_data, 
    radius=1000  # 1 km minimum separation
)
print(f"Removed {len(clean_data) - len(final_data)} duplicate stations")

# Now interpolate with cleaned data
grid_x, grid_y, clean_grid = mpinterp.natural_neighbor_to_grid(
    final_x, final_y, final_data,
    hres=10000,
    minimum_neighbors=2
)

# Atmospheric sounding interpolation to standard levels
pressure_obs = np.array([1000, 925, 850, 700, 500, 400, 300, 250]) * units.hPa
temp_obs = np.array([20, 18, 15, 8, -5, -15, -30, -40]) * units.celsius

# Standard pressure levels
standard_levels = np.array([1000, 925, 850, 700, 500, 400, 300, 250, 200, 150, 100]) * units.hPa

# Interpolate to standard levels
temp_interp = mpinterp.interpolate_1d(
    standard_levels.m, pressure_obs.m, temp_obs.m,
    fill_value=np.nan  # Don't extrapolate beyond observed range
)

print("Original levels:", len(pressure_obs))
print("Standard levels:", len(standard_levels))
print("Interpolated temperatures:", temp_interp)

# Plot sounding
plt.figure(figsize=(8, 10))
plt.semilogy(temp_obs, pressure_obs, 'ro-', label='Observed', markersize=8)
plt.semilogy(temp_interp, standard_levels, 'b*-', label='Interpolated', markersize=6)
plt.gca().invert_yaxis()
plt.xlabel('Temperature (°C)')
plt.ylabel('Pressure (hPa)')
plt.title('Atmospheric Sounding - Standard Level Interpolation')
plt.legend()
plt.grid(True)
plt.show()

import xarray as xr
import metpy.xarray

# Load 3D atmospheric data
ds = xr.open_dataset('gfs_analysis.nc').metpy.parse_cf()
temperature_3d = ds['temperature']

# Define cross-section path
start_point = [35.0, -100.0]  # [lat, lon] 
end_point = [45.0, -90.0]     # [lat, lon]

# Extract cross-section
cross_sec = mpinterp.cross_section(
    temperature_3d, 
    start=start_point, 
    end=end_point,
    steps=50,
    interp_type='linear'
)

# Plot cross-section
plt.figure(figsize=(12, 8))
plt.contourf(cross_sec.distance, cross_sec.pressure, cross_sec.T, 
            levels=20, cmap='RdYlBu_r')
plt.colorbar(label='Temperature (K)')
plt.gca().invert_yaxis()
plt.xlabel('Distance along cross-section')
plt.ylabel('Pressure (hPa)')
plt.title('Temperature Cross-Section')
plt.show()

# Working with lat/lon station data
station_lons = np.array([-105, -100, -95, -90, -85]) * units.degrees_east
station_lats = np.array([40, 41, 39, 42, 38]) * units.degrees_north
dewpoints = np.array([10, 8, 12, 6, 14]) * units.celsius

# Define analysis domain 
boundary = {
    'west': -110 * units.degrees_east,
    'east': -80 * units.degrees_east, 
    'south': 35 * units.degrees_north,
    'north': 45 * units.degrees_north
}

# Interpolate from geographic to projected grid
from metpy.crs import LambertConformal
proj = LambertConformal(central_longitude=-95, central_latitude=40)

grid_x, grid_y, dewpoint_grid = mpinterp.geodetic_to_grid(
    station_lons.m, station_lats.m, dewpoints.m,
    boundary_coords=boundary,
    hres=50000,  # 50 km in projected coordinates
    grid_projection=proj,
    interp_type='linear'
)

print(f"Grid shape: {dewpoint_grid.shape}")
print(f"Longitude range: {station_lons.min()} to {station_lons.max()}")
print(f"Latitude range: {station_lats.min()} to {station_lats.max()}")

# Time series with gaps
time_series = np.array([15.0, 16.2, np.nan, np.nan, 18.5, 19.1, np.nan, 20.3])
times = np.arange(len(time_series))

# Fill missing values with interpolation
filled_series = mpinterp.interpolate_nans(
    time_series, 
    method='linear',
    limit=2  # Only fill gaps of 2 or fewer points
)

# Plot original and filled data
plt.figure(figsize=(10, 6))
plt.plot(times, time_series, 'ro-', label='Original (with NaNs)', markersize=8)
plt.plot(times, filled_series, 'b*-', label='Interpolated', markersize=6)
plt.xlabel('Time Index')
plt.ylabel('Temperature (°C)')
plt.title('Gap Filling with Linear Interpolation')
plt.legend()
plt.grid(True)
plt.show()

print("Original:", time_series)
print("Filled:  ", filled_series)

from metpy.io import parse_metar_to_dataframe
import metpy.interpolate as mpinterp

# Load surface observations
df = parse_metar_to_dataframe('surface_obs.txt')

# Extract coordinates and temperature
lons = df['longitude'].values
lats = df['latitude'].values  
temps = df['air_temperature'].values

# Remove missing data
clean_lons, clean_lats, clean_temps = mpinterp.remove_nan_observations(
    lons, lats, temps
)

# Create temperature analysis
grid_x, grid_y, temp_analysis = mpinterp.natural_neighbor_to_grid(
    clean_lons, clean_lats, clean_temps,
    hres=50000,
    interp_type='linear'
)

import metpy.calc as mpcalc

# After interpolation, use with MetPy calculations
# Convert grid back to DataArray with coordinates
import xarray as xr

temp_da = xr.DataArray(
    temp_analysis * units.celsius,
    dims=['y', 'x'],
    coords={'y': grid_y, 'x': grid_x}
)

# Now use with MetPy calculations
# Calculate temperature gradient
temp_grad = mpcalc.gradient(temp_da)
print("Temperature gradient calculated from interpolated field")

from typing import Tuple, Optional, Union, Sequence
import numpy as np
import xarray as xr
from pint import Quantity

# Input data types
CoordinateArray = Union[np.ndarray, Sequence, Quantity]
DataArray = Union[np.ndarray, Sequence, Quantity]

# Grid output types
GridTuple = Tuple[np.ndarray, np.ndarray, np.ndarray]  # (x, y, data)
CrossSection = xr.DataArray

# Parameter types
InterpolationType = str  # 'linear', 'nearest', 'cubic'
SearchRadius = Union[float, int, Quantity]
Resolution = Union[float, int, Quantity]

# Boundary specification
BoundaryCoords = Union[Dict[str, float], Sequence[float]]