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