1# Copyright (c) 2016 MetPy Developers.
2# Distributed under the terms of the BSD 3-Clause License.
3# SPDX-License-Identifier: BSD-3-Clause
4"""
5Point Interpolation
6===================
7
8Compares different point interpolation approaches.
9"""
10
11import cartopy.crs as ccrs
12import cartopy.feature as cfeature
13from matplotlib.colors import BoundaryNorm
14import matplotlib.pyplot as plt
15import numpy as np
16
17from metpy.cbook import get_test_data
18from metpy.interpolate import (interpolate_to_grid, remove_nan_observations,
19                               remove_repeat_coordinates)
20from metpy.plots import add_metpy_logo
21
22
23###########################################
24def basic_map(proj, title):
25    """Make our basic default map for plotting"""
26    fig = plt.figure(figsize=(15, 10))
27    add_metpy_logo(fig, 0, 80, size='large')
28    view = fig.add_axes([0, 0, 1, 1], projection=proj)
29    view.set_title(title)
30    view.set_extent([-120, -70, 20, 50])
31    view.add_feature(cfeature.STATES.with_scale('50m'))
32    view.add_feature(cfeature.OCEAN)
33    view.add_feature(cfeature.COASTLINE)
34    view.add_feature(cfeature.BORDERS, linestyle=':')
35    return fig, view
36
37
38def station_test_data(variable_names, proj_from=None, proj_to=None):
39    with get_test_data('station_data.txt') as f:
40        all_data = np.loadtxt(f, skiprows=1, delimiter=',',
41                              usecols=(1, 2, 3, 4, 5, 6, 7, 17, 18, 19),
42                              dtype=np.dtype([('stid', '3S'), ('lat', 'f'), ('lon', 'f'),
43                                              ('slp', 'f'), ('air_temperature', 'f'),
44                                              ('cloud_fraction', 'f'), ('dewpoint', 'f'),
45                                              ('weather', '16S'),
46                                              ('wind_dir', 'f'), ('wind_speed', 'f')]))
47
48    all_stids = [s.decode('ascii') for s in all_data['stid']]
49
50    data = np.concatenate([all_data[all_stids.index(site)].reshape(1, ) for site in all_stids])
51
52    value = data[variable_names]
53    lon = data['lon']
54    lat = data['lat']
55
56    if proj_from is not None and proj_to is not None:
57
58        try:
59
60            proj_points = proj_to.transform_points(proj_from, lon, lat)
61            return proj_points[:, 0], proj_points[:, 1], value
62
63        except Exception as e:
64
65            print(e)
66            return None
67
68    return lon, lat, value
69
70
71from_proj = ccrs.Geodetic()
72to_proj = ccrs.AlbersEqualArea(central_longitude=-97.0000, central_latitude=38.0000)
73
74levels = list(range(-20, 20, 1))
75cmap = plt.get_cmap('magma')
76norm = BoundaryNorm(levels, ncolors=cmap.N, clip=True)
77
78x, y, temp = station_test_data('air_temperature', from_proj, to_proj)
79
80x, y, temp = remove_nan_observations(x, y, temp)
81x, y, temp = remove_repeat_coordinates(x, y, temp)
82
83###########################################
84# Scipy.interpolate linear
85# ------------------------
86gx, gy, img = interpolate_to_grid(x, y, temp, interp_type='linear', hres=75000)
87img = np.ma.masked_where(np.isnan(img), img)
88fig, view = basic_map(to_proj, 'Linear')
89mmb = view.pcolormesh(gx, gy, img, cmap=cmap, norm=norm)
90fig.colorbar(mmb, shrink=.4, pad=0, boundaries=levels)
91
92###########################################
93# Natural neighbor interpolation (MetPy implementation)
94# -----------------------------------------------------
95# `Reference <https://github.com/Unidata/MetPy/files/138653/cwp-657.pdf>`_
96gx, gy, img = interpolate_to_grid(x, y, temp, interp_type='natural_neighbor', hres=75000)
97img = np.ma.masked_where(np.isnan(img), img)
98fig, view = basic_map(to_proj, 'Natural Neighbor')
99mmb = view.pcolormesh(gx, gy, img, cmap=cmap, norm=norm)
100fig.colorbar(mmb, shrink=.4, pad=0, boundaries=levels)
101
102###########################################
103# Cressman interpolation
104# ----------------------
105# search_radius = 100 km
106#
107# grid resolution = 25 km
108#
109# min_neighbors = 1
110gx, gy, img = interpolate_to_grid(x, y, temp, interp_type='cressman', minimum_neighbors=1,
111                                  hres=75000, search_radius=100000)
112img = np.ma.masked_where(np.isnan(img), img)
113fig, view = basic_map(to_proj, 'Cressman')
114mmb = view.pcolormesh(gx, gy, img, cmap=cmap, norm=norm)
115fig.colorbar(mmb, shrink=.4, pad=0, boundaries=levels)
116
117###########################################
118# Barnes Interpolation
119# --------------------
120# search_radius = 100km
121#
122# min_neighbors = 3
123gx, gy, img1 = interpolate_to_grid(x, y, temp, interp_type='barnes', hres=75000,
124                                   search_radius=100000)
125img1 = np.ma.masked_where(np.isnan(img1), img1)
126fig, view = basic_map(to_proj, 'Barnes')
127mmb = view.pcolormesh(gx, gy, img1, cmap=cmap, norm=norm)
128fig.colorbar(mmb, shrink=.4, pad=0, boundaries=levels)
129
130###########################################
131# Radial basis function interpolation
132# ------------------------------------
133# linear
134gx, gy, img = interpolate_to_grid(x, y, temp, interp_type='rbf', hres=75000, rbf_func='linear',
135                                  rbf_smooth=0)
136img = np.ma.masked_where(np.isnan(img), img)
137fig, view = basic_map(to_proj, 'Radial Basis Function')
138mmb = view.pcolormesh(gx, gy, img, cmap=cmap, norm=norm)
139fig.colorbar(mmb, shrink=.4, pad=0, boundaries=levels)
140
141plt.show()
142