0
# Spatial Interpolation
1

2
MetPy provides advanced interpolation functions specifically designed for meteorological data analysis. These tools handle sparse station observations, irregular grids, and three-dimensional atmospheric data with sophisticated methods including natural neighbor interpolation, objective analysis techniques, and cross-section extraction.
3

4
## Capabilities
5

6
### Grid-Based Interpolation
7

8
Transform sparse station observations to regular grids using various interpolation methods optimized for meteorological applications.
9

10
```python { .api }
11
def natural_neighbor_to_grid(xp, yp, variable, interp_type='linear', hres=50000, 
12
                           search_radius=None, rbf_func='linear', rbf_smooth=0,
13
                           minimum_neighbors=3, gamma=0.25, kappa_star=5.052, 
14
                           boundary_coords=None, grid_projection=None):
15
    """
16
    Generate a natural neighbor interpolation of the given point data to a grid.
17
    
18
    Natural neighbor interpolation is a method of spatial interpolation for
19
    multivariate data that preserves local properties and is well-suited for
20
    irregular station networks common in meteorology.
21
    
22
    Parameters:
23
    - xp: x-coordinates of observation points
24
    - yp: y-coordinates of observation points  
25
    - variable: data values at observation points
26
    - interp_type: interpolation method ('linear', 'nearest', 'cubic')
27
    - hres: horizontal resolution of output grid
28
    - search_radius: maximum search radius for points
29
    - rbf_func: radial basis function type
30
    - rbf_smooth: smoothing parameter for RBF
31
    - minimum_neighbors: minimum neighbors required for interpolation
32
    - gamma: shape parameter for natural neighbor weighting
33
    - kappa_star: scaling parameter for influence radius
34
    - boundary_coords: boundary coordinates for interpolation domain
35
    - grid_projection: coordinate reference system for output grid
36
    
37
    Returns:
38
    Tuple of (grid_x, grid_y, grid_data) arrays with interpolated values
39
    """
40

41
def inverse_distance_to_grid(xp, yp, variable, hres, search_radius=None, 
42
                           rbf_func='linear', rbf_smooth=0, minimum_neighbors=3,
43
                           gamma=0.25, kappa_star=5.052, boundary_coords=None):
44
    """
45
    Interpolate point data to grid using inverse distance weighting.
46
    
47
    Inverse distance weighting is a simple but effective interpolation method
48
    that weights observations inversely proportional to their distance from
49
    the interpolation point.
50
    
51
    Parameters:
52
    - xp: x-coordinates of observation points
53
    - yp: y-coordinates of observation points
54
    - variable: data values at observation points
55
    - hres: horizontal resolution of output grid
56
    - search_radius: maximum search radius for points
57
    - rbf_func: radial basis function type
58
    - rbf_smooth: smoothing parameter
59
    - minimum_neighbors: minimum neighbors for interpolation
60
    - gamma: weighting exponent (typically 1-3)
61
    - kappa_star: influence radius scaling
62
    - boundary_coords: interpolation domain boundary
63
    
64
    Returns:
65
    Tuple of (grid_x, grid_y, grid_data) arrays
66
    """
67

68
def barnes_to_grid(xp, yp, variable, hres, search_radius, kappa_star=5.052,
69
                  gamma=0.25, minimum_neighbors=3, boundary_coords=None):
70
    """
71
    Interpolate point data using Barnes objective analysis.
72
    
73
    Barnes analysis is a successive correction method widely used in meteorology
74
    for creating gridded analyses from irregular station data. It applies
75
    Gaussian weighting with multiple passes to reduce analysis error.
76
    
77
    Parameters:
78
    - xp: x-coordinates of observation points
79
    - yp: y-coordinates of observation points
80
    - variable: data values at observation points
81
    - hres: horizontal resolution of output grid
82
    - search_radius: search radius for station selection
83
    - kappa_star: Gaussian response parameter (typically 3-10)
84
    - gamma: convergence parameter for successive corrections
85
    - minimum_neighbors: minimum stations required for analysis
86
    - boundary_coords: analysis domain boundary coordinates
87
    
88
    Returns:
89
    Tuple of (grid_x, grid_y, grid_data) with Barnes analysis
90
    """
91

92
def cressman_to_grid(xp, yp, variable, hres, search_radius, minimum_neighbors=3,
93
                    boundary_coords=None):
94
    """
95
    Interpolate point data using Cressman objective analysis.
96
    
97
    Cressman analysis uses successive corrections with distance-weighted
98
    influence functions. It's computationally efficient and widely used
99
    in operational meteorology for real-time analysis.
100
    
101
    Parameters:
102
    - xp: x-coordinates of observation points
103
    - yp: y-coordinates of observation points
104
    - variable: data values at observation points
105
    - hres: horizontal resolution of output grid
106
    - search_radius: influence radius for each correction pass
107
    - minimum_neighbors: minimum observations for grid point analysis
108
    - boundary_coords: analysis domain boundaries
109
    
110
    Returns:
111
    Tuple of (grid_x, grid_y, grid_data) with Cressman analysis
112
    """
113

114
def geodetic_to_grid(xp, yp, variable, boundary_coords, hres, interp_type='linear',
115
                    search_radius=None, rbf_func='linear', rbf_smooth=0,
116
                    minimum_neighbors=3, gamma=0.25, grid_projection=None):
117
    """
118
    Interpolate data from geodetic coordinates (lat/lon) to a projected grid.
119
    
120
    Handles coordinate transformation from geographic to projected coordinate
121
    systems while performing spatial interpolation, essential for meteorological
122
    analysis spanning large geographic areas.
123
    
124
    Parameters:
125
    - xp: longitude coordinates of observation points
126
    - yp: latitude coordinates of observation points
127
    - variable: data values at observation points
128
    - boundary_coords: geographic bounds for interpolation domain
129
    - hres: horizontal resolution in projected coordinates
130
    - interp_type: interpolation method
131
    - search_radius: search radius in projected coordinates
132
    - rbf_func: radial basis function type
133
    - rbf_smooth: smoothing parameter
134
    - minimum_neighbors: minimum neighbors for interpolation
135
    - gamma: method-specific parameter
136
    - grid_projection: target projection for output grid
137
    
138
    Returns:
139
    Tuple of (grid_x, grid_y, grid_data) in projected coordinates
140
    """
141
```
142

143
### Data Quality and Preprocessing
144

145
Utilities for cleaning and preparing observational data for interpolation.
146

147
```python { .api }
148
def remove_nan_observations(x, y, *variables):
149
    """
150
    Remove observations with NaN values from coordinate and data arrays.
151
    
152
    Quality control step essential for reliable interpolation results,
153
    removing invalid or missing observations that would contaminate
154
    the analysis.
155
    
156
    Parameters:
157
    - x: x-coordinate array
158
    - y: y-coordinate array  
159
    - variables: one or more data variable arrays
160
    
161
    Returns:
162
    Tuple of cleaned (x, y, *variables) with NaN observations removed
163
    """
164

165
def remove_repeat_coordinates(x, y, z, radius=None):
166
    """
167
    Remove observations with duplicate or nearly duplicate coordinates.
168
    
169
    Eliminates station observations that are too close together, which
170
    can cause numerical issues in interpolation algorithms and bias
171
    the analysis toward over-sampled regions.
172
    
173
    Parameters:
174
    - x: x-coordinate array
175
    - y: y-coordinate array
176
    - z: data value array
177
    - radius: minimum separation distance between stations
178
    
179
    Returns:
180
    Tuple of (x, y, z) with duplicate coordinates removed
181
    """
182
```
183

184
### One-Dimensional Interpolation
185

186
Linear interpolation functions for vertical profiles and time series.
187

188
```python { .api }
189
def interpolate_1d(x, xp, *variables, axis=0, fill_value=np.nan):
190
    """
191
    Interpolate 1D arrays to new coordinate values.
192
    
193
    Essential for vertical interpolation of atmospheric soundings to
194
    standard pressure levels, height coordinates, or potential temperature
195
    surfaces commonly used in meteorological analysis.
196
    
197
    Parameters:
198
    - x: new coordinate values for interpolation
199
    - xp: original coordinate values (must be monotonic)
200
    - variables: one or more data arrays to interpolate
201
    - axis: axis along which to interpolate
202
    - fill_value: value for extrapolation outside data range
203
    
204
    Returns:
205
    Interpolated arrays at new coordinate values
206
    """
207

208
def interpolate_nans(array, method='linear', axis=0, limit=None):
209
    """
210
    Fill NaN values in array using interpolation.
211
    
212
    Useful for filling missing values in time series or vertical profiles
213
    while preserving the overall structure and variability of the data.
214
    
215
    Parameters:
216
    - array: input array with NaN values to fill
217
    - method: interpolation method ('linear', 'nearest', 'cubic')
218
    - axis: axis along which to interpolate
219
    - limit: maximum number of consecutive NaNs to fill
220
    
221
    Returns:
222
    Array with NaN values filled by interpolation
223
    """
224
```
225

226
### Cross-Section Analysis
227

228
Extract and interpolate data along arbitrary paths through three-dimensional atmospheric data.
229

230
```python { .api }
231
def cross_section(data, start, end, steps=100, interp_type='linear'):
232
    """
233
    Extract a cross-section through gridded 3D atmospheric data.
234
    
235
    Creates vertical cross-sections along arbitrary horizontal paths,
236
    essential for analyzing atmospheric structure, frontal zones,
237
    and three-dimensional meteorological phenomena.
238
    
239
    Parameters:
240
    - data: 3D xarray DataArray with atmospheric data
241
    - start: starting point coordinates (lat, lon) or (x, y)
242
    - end: ending point coordinates (lat, lon) or (x, y)
243
    - steps: number of points along cross-section path
244
    - interp_type: interpolation method for extracting values
245
    
246
    Returns:
247
    xarray DataArray containing the cross-section with distance and
248
    vertical coordinates
249
    """
250
```
251

252
## Usage Examples
253

254
### Basic Grid Interpolation
255

256
```python
257
import metpy.interpolate as mpinterp
258
from metpy.units import units
259
import numpy as np
260
import matplotlib.pyplot as plt
261

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

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

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

278
# Plot results
279
plt.figure(figsize=(10, 8))
280
plt.contourf(grid_x, grid_y, temp_grid, levels=15, cmap='RdYlBu_r')
281
plt.colorbar(label='Temperature (°C)')
282
plt.scatter(x_coords, y_coords, c='black', s=50, label='Stations')
283
plt.xlabel('X (m)')
284
plt.ylabel('Y (m)')
285
plt.title('Temperature Analysis - Natural Neighbor Interpolation')
286
plt.legend()
287
plt.show()
288
```
289

290
### Barnes Objective Analysis
291

292
```python
293
# More sophisticated analysis with Barnes method
294
search_radius = 200000  # 200 km search radius
295

296
grid_x, grid_y, temp_barnes = mpinterp.barnes_to_grid(
297
    x_coords, y_coords, temperatures.m,
298
    hres=25000,  # 25 km resolution
299
    search_radius=search_radius,
300
    kappa_star=5.052,  # Standard Barnes parameter
301
    gamma=0.25,        # Convergence parameter
302
    minimum_neighbors=3
303
)
304

305
# Compare with Cressman analysis
306
grid_x_cress, grid_y_cress, temp_cressman = mpinterp.cressman_to_grid(
307
    x_coords, y_coords, temperatures.m,
308
    hres=25000,
309
    search_radius=search_radius,
310
    minimum_neighbors=3
311
)
312

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

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

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

328
plt.tight_layout()
329
plt.show()
330
```
331

332
### Data Quality Control
333

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

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

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

353
# Now interpolate with cleaned data
354
grid_x, grid_y, clean_grid = mpinterp.natural_neighbor_to_grid(
355
    final_x, final_y, final_data,
356
    hres=10000,
357
    minimum_neighbors=2
358
)
359
```
360

361
### Vertical Profile Interpolation
362

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

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

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

377
print("Original levels:", len(pressure_obs))
378
print("Standard levels:", len(standard_levels))
379
print("Interpolated temperatures:", temp_interp)
380

381
# Plot sounding
382
plt.figure(figsize=(8, 10))
383
plt.semilogy(temp_obs, pressure_obs, 'ro-', label='Observed', markersize=8)
384
plt.semilogy(temp_interp, standard_levels, 'b*-', label='Interpolated', markersize=6)
385
plt.gca().invert_yaxis()
386
plt.xlabel('Temperature (°C)')
387
plt.ylabel('Pressure (hPa)')
388
plt.title('Atmospheric Sounding - Standard Level Interpolation')
389
plt.legend()
390
plt.grid(True)
391
plt.show()
392
```
393

394
### Cross-Section Analysis
395

396
```python
397
import xarray as xr
398
import metpy.xarray
399

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

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

408
# Extract cross-section
409
cross_sec = mpinterp.cross_section(
410
    temperature_3d, 
411
    start=start_point, 
412
    end=end_point,
413
    steps=50,
414
    interp_type='linear'
415
)
416

417
# Plot cross-section
418
plt.figure(figsize=(12, 8))
419
plt.contourf(cross_sec.distance, cross_sec.pressure, cross_sec.T, 
420
            levels=20, cmap='RdYlBu_r')
421
plt.colorbar(label='Temperature (K)')
422
plt.gca().invert_yaxis()
423
plt.xlabel('Distance along cross-section')
424
plt.ylabel('Pressure (hPa)')
425
plt.title('Temperature Cross-Section')
426
plt.show()
427
```
428

429
### Geographic Coordinate Interpolation
430

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

437
# Define analysis domain 
438
boundary = {
439
    'west': -110 * units.degrees_east,
440
    'east': -80 * units.degrees_east, 
441
    'south': 35 * units.degrees_north,
442
    'north': 45 * units.degrees_north
443
}
444

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

449
grid_x, grid_y, dewpoint_grid = mpinterp.geodetic_to_grid(
450
    station_lons.m, station_lats.m, dewpoints.m,
451
    boundary_coords=boundary,
452
    hres=50000,  # 50 km in projected coordinates
453
    grid_projection=proj,
454
    interp_type='linear'
455
)
456

457
print(f"Grid shape: {dewpoint_grid.shape}")
458
print(f"Longitude range: {station_lons.min()} to {station_lons.max()}")
459
print(f"Latitude range: {station_lats.min()} to {station_lats.max()}")
460
```
461

462
### Handling Missing Data
463

464
```python
465
# Time series with gaps
466
time_series = np.array([15.0, 16.2, np.nan, np.nan, 18.5, 19.1, np.nan, 20.3])
467
times = np.arange(len(time_series))
468

469
# Fill missing values with interpolation
470
filled_series = mpinterp.interpolate_nans(
471
    time_series, 
472
    method='linear',
473
    limit=2  # Only fill gaps of 2 or fewer points
474
)
475

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

487
print("Original:", time_series)
488
print("Filled:  ", filled_series)
489
```
490

491
## Integration with MetPy Ecosystem
492

493
### Working with MetPy I/O
494

495
```python
496
from metpy.io import parse_metar_to_dataframe
497
import metpy.interpolate as mpinterp
498

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

502
# Extract coordinates and temperature
503
lons = df['longitude'].values
504
lats = df['latitude'].values  
505
temps = df['air_temperature'].values
506

507
# Remove missing data
508
clean_lons, clean_lats, clean_temps = mpinterp.remove_nan_observations(
509
    lons, lats, temps
510
)
511

512
# Create temperature analysis
513
grid_x, grid_y, temp_analysis = mpinterp.natural_neighbor_to_grid(
514
    clean_lons, clean_lats, clean_temps,
515
    hres=50000,
516
    interp_type='linear'
517
)
518
```
519

520
### Integration with Calculations
521

522
```python
523
import metpy.calc as mpcalc
524

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

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

535
# Now use with MetPy calculations
536
# Calculate temperature gradient
537
temp_grad = mpcalc.gradient(temp_da)
538
print("Temperature gradient calculated from interpolated field")
539
```
540

541
## Types
542

543
```python { .api }
544
from typing import Tuple, Optional, Union, Sequence
545
import numpy as np
546
import xarray as xr
547
from pint import Quantity
548

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

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

557
# Parameter types
558
InterpolationType = str  # 'linear', 'nearest', 'cubic'
559
SearchRadius = Union[float, int, Quantity]
560
Resolution = Union[float, int, Quantity]
561

562
# Boundary specification
563
BoundaryCoords = Union[Dict[str, float], Sequence[float]]
564
```