title: Point Interpolation with Metpython and Basemap
body: <p>I change the <a href="https://unidata.github.io/MetPy/latest/examples/gridding/Point_Interpolation.html" rel="nofollow noreferrer">script</a>, and try to show it by using basemap.</p>
  <pre><code># Copyright (c) 2016 MetPy Developers.
  # Distributed under the terms of the BSD 3-Clause License.
  # SPDX-License-Identifier: BSD-3-Clause
  &quot;&quot;&quot;   
  [Point_Interpolation](https://unidata.github.io/MetPy/latest/examples/gridding/Point_Interpolation.html?highlight=basic_map)
  
  &quot;&quot;&quot;
  </code></pre>
  <p>Import some library:</p>
  <pre><code>import cartopy.crs as ccrs
  import cartopy.feature as cfeature
  from matplotlib.colors import BoundaryNorm
  import matplotlib.pyplot as plt
  from matplotlib import cm as CM
  from mpl_toolkits.basemap import Basemap
  from matplotlib.patches import Polygon
  from metpy.cbook import get_test_data
  from metpy.gridding.gridding_functions import (interpolate, remove_nan_observations,
                                                 remove_repeat_coordinates)
  import pandas as pd
  import numpy as np
  import numpy.ma as ma
  </code></pre>
  <p>Here some fuction:</p>
  <pre><code>def basic_map(proj):
      &quot;&quot;&quot;Make our basic default map for plotting&quot;&quot;&quot;
      fig = plt.figure(figsize=(15, 10))
      #add_metpy_logo(fig, 0, 80, size='large')
      view = fig.add_axes([0, 0, 1, 1], projection=proj)
      view.set_extent([100, 130, 20, 50])
      view.add_feature(cfeature.STATES.with_scale('50m'))
      view.add_feature(cfeature.OCEAN)
      view.add_feature(cfeature.COASTLINE)
      view.add_feature(cfeature.BORDERS, linestyle=':')
      return fig, view
  </code></pre>
  <p>get station_test_data:</p>
  <pre><code>def station_test_data(variable_names, proj_from=None, proj_to=None):
      #with get_test_data('station_data.txt') as f:
      all_data = np.loadtxt('2016072713.000', skiprows=1,
      #     all_data = np.loadtxt(f, skiprows=1,delimiter=',',
      
                                usecols=(0,1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11),
                                dtype=np.dtype([('stid', '5S'), ('lat', 'f'), ('lon', 'f'),
                                                ('aqi', 'f'), ('grade', 'f'),
                                                ('pm25', 'f'), ('pm10', 'f'),
                                                ('co', 'f'),('no2','f'),('o3','f'),
                                                ('o38h', 'f'), ('so2', 'f')]))
     
      all_stids = [s.decode('ascii') for s in all_data['stid']]
      data = np.concatenate([all_data[all_stids.index(site)].reshape(1, ) for site in all_stids])
  
      value = data[variable_names]
      lon = data['lon']
      lat = data['lat']
  
      if proj_from is not None and proj_to is not None:
  
          try:
  
              proj_points = proj_to.transform_points(proj_from, lon, lat)
              return proj_points[:, 0], proj_points[:, 1], value
  
          except Exception as e:
  
              print(e)
              return None
  
      return lon, lat, value
  def remove_inf_x(x, y, z):
      x_ = x[~np.isinf(x)]
      y_ = y[~np.isinf(x)]
      z_ = z[~np.isinf(x)]
  
      return x_, y_, z_
  def remove_inf_y(x, y, z):
      x_ = x[~np.isinf(y)]
      y_ = y[~np.isinf(y)]
      z_ = z[~np.isinf(y)]
  
      return x_, y_, z_
  </code></pre>
  <p>add ploygen plot</p>
  <pre><code>def plot_rec(map,lower_left,upper_left,upper_right,lower_right,text):
      lon_poly = np.array([lower_left[0], upper_left[0],upper_right[0], lower_right[0]])
      lat_poly = np.array([lower_left[1], upper_left[1],upper_right[1], lower_right[1]])
      X, Y  = map(lon_poly, lat_poly)
      xy = np.vstack([X,Y]).T
  
      poly = Polygon(xy, closed=True, \
              facecolor='None',edgecolor='red', \
              linewidth=1.\
              )
      ax, ay = map(lower_left[0],lower_left[1])
      plt.text(ax, ay,text,fontsize=6,fontweight='bold',
                      ha='left',va='bottom',color='k')
      plt.gca().add_patch(poly)
  </code></pre>
  <p>Get the data:</p>
  <pre><code>from_proj = ccrs.Geodetic()
  to_proj = ccrs.AlbersEqualArea(central_longitude=110.0000, central_latitude=32.0000)
  x, y, temp = station_test_data('pm25', from_proj, to_proj)
  
  x, y, temp = remove_nan_observations(x, y, temp)
  x, y, temp = remove_inf_x(x, y, temp)
  x, y, temp = remove_inf_y(x, y, temp)
  x, y, temp = remove_repeat_coordinates(x, y, temp)
  
  ###########################################
  </code></pre>
  <p>Barnes Interpolation: search_radius = 100km min_neighbors = 3</p>
  <pre><code>gx, gy, img1 = interpolate(x, y, temp, interp_type='barnes', hres=75000, search_radius=100000)
  img1 = np.ma.masked_where(np.isnan(img1), img1)
  </code></pre>
  <p>Then I make the plot:</p>
  <pre><code>fig = plt.figure(figsize=(8, 8))
  #fig, view = basic_map(to_proj)
  m = Basemap(width=8000000,height=7000000,
              resolution='l',projection='aea',\
              lat_1=0.,lat_2=40,lon_0=110,lat_0=20)
  #m.shadedrelief()
  m.drawparallels(np.arange(20.,40,2.5),linewidth=1, dashes=[4, 2], labels=[1,0,0,0], color= 'gray',zorder=0, fontsize=10)
  m.drawmeridians(np.arange(100.,125.,2.),linewidth=1, dashes=[4, 2], labels=[0,0,0,1], color= 'gray',zorder=0, fontsize=10)
  m.readshapefile('dijishi_2004','state',color='gray')
  levels = list(range(0, 200, 10))
  cmap = plt.get_cmap('viridis')
  norm = BoundaryNorm(levels, ncolors=cmap.N, clip=True)
  mmb = m.pcolormesh(gx, gy, img1, cmap=cmap, norm=norm)
  plt.colorbar(label=r'$PM_{2.5}(\mu g/m^3 $)') 
  plt.show()
  </code></pre>
  <p><a href="https://i.stack.imgur.com/AnO7x.jpg" rel="nofollow noreferrer">It seems the location is wrong!</a></p>
  <p>Maybe I use the Point_Interpolation, which has changed the coordinate.</p>
  <p>Here is the import data which I use like this:</p>
  <blockquote>
  <p>99306 31.1654 121.412 NaN NaN NaN 37.000 0.875 10.000 141.000 NaN NaN</p>
  <p>99299 31.2036 121.478 NaN NaN NaN 49.000 0.420 18.000 157.000 NaN 16.000</p>
  <p>99302 31.1907 121.703 NaN NaN 75.000 112.000 1.571 54.000 167.000 NaN 34.000</p>
  <p>99300 31.3008 121.467 NaN NaN 53.000 NaN 0.414 21.000 128.000 NaN 10.000</p>
  <p>99304 31.2071 121.577 NaN NaN 20.000 66.000 NaN 20.000 192.000 NaN 28.000</p>
  <p>99305 31.0935 120.978 NaN NaN NaN 5.000 0.717 23.000 140.000 NaN 13.000</p>
  </blockquote>