-1

我目前正在尝试使用本网站提供的代码（https://unidata.github.io/MetPy/latest/examples/gridding/Point_Interpolation.html#sphx-glr-examples-gridding-point-interpolation-py）在 Jupyter 笔记本上使用数据的线性插值创建台湾地图。

我的数据是这种形式：

17070123, lat, lon, tem

C0A92, 25.27, 121.56, 29.3

C0AD0, 25.26, 121.49, 28.2

C0A94, 25.23, 121.64, 26.2

46691, 25.19, 121.52, 23.4

46690, 25.17, 121.44, 27.3

46693, 25.17, 121.54, 22.5

C0AD1, 25.15, 121.4, 28.5

46694, 25.13, 121.73, 28.6

C0A95, 25.13, 121.92, -999

C0A9B, 25.12, 121.51, 26.8

C0A9C, 25.12, 121.53, 28.3

C0A66, 25.11, 121.79, 27.8

C0A98, 25.11, 121.46, 29.6

C0A68, 25.09, 121.43, -999

也以这种形式：

#17070123   lat lon T

C0A92   25.27   121.56  29.3

C0AD0   25.26   121.49  28.2

C0A94   25.23   121.64  26.2

46691   25.19   121.52  23.4

46690   25.17   121.44  27.3

46693   25.17   121.54  22.5

C0AD1   25.15   121.4   28.5

46694   25.13   121.73  28.6

C0A95   25.13   121.92  -999

C0A9B   25.12   121.51  26.8

C0A9C   25.12   121.53  28.3

C0A66   25.11   121.79  27.8

C0A98   25.11   121.46  29.6

C0A68   25.09   121.43  -999

我的代码如下所示：

# In[1]:

import cartopy
import cartopy.crs as ccrs
from matplotlib.colors import BoundaryNorm
import matplotlib.pyplot as plt
import numpy as np


# In[2]:

from metpy.cbook import get_test_data
from metpy.gridding.gridding_functions import (interpolate, 
remove_nan_observations,
                                           remove_repeat_coordinates)


# In[3]:

def basic_map(map_proj):
   """Make our basic default map for plotting"""
   fig = plt.figure(figsize=(15, 10))
   view = fig.add_axes([0, 0, 1, 1], projection=to_proj)
   view.set_extent([120.5, 122.5, 24.5, 25.5])
   view.add_feature(cartopy.feature.NaturalEarthFeature(category='cultural',

   name='admin_1_states_provinces_lakes',
                                                    scale='50m', 
facecolor='none'))
   view.add_feature(cartopy.feature.OCEAN)
   view.add_feature(cartopy.feature.COASTLINE)
   view.add_feature(cartopy.feature.BORDERS, linestyle=':')
   return view


# In[4]:

def station_test_data(variable_names, proj_from=None, proj_to=None):
    f = ('temp.txt')
    all_data = np.loadtxt(f, skiprows=0, delimiter='\t',
                      usecols=(0, 1, 2, 3),
                      dtype=np.dtype([('stid', '5S'), ('lat', 'f'), ('lon', 
                       'f'), ('T', '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


# In[5]:

from_proj = ccrs.Geodetic()
to_proj = ccrs.AlbersEqualArea(central_longitude=120.0000, 
central_latitude=25.0000)


# In[6]:

levels = list(range(20, 30, 1))
cmap = plt.get_cmap('magma')
norm = BoundaryNorm(levels, ncolors=cmap.N, clip=True)


# In[7]:

x, y, temp = station_test_data('T', from_proj, to_proj)


# In[8]:

x, y, temp = remove_nan_observations(x, y, temp)
x, y, temp = remove_repeat_coordinates(x, y, temp)


# In[9]:

gx, gy, img = interpolate(x, y, temp, interp_type='linear', hres=75000)
img = np.ma.masked_where(np.isnan(img), img)
view = basic_map(to_proj)
mmb = view.pcolormesh(gx, gy, img, cmap=cmap, norm=norm)
plt.colorbar(mmb, shrink=.4, pad=0, boundaries=levels)


# In[10]:

#Show map of TW with interpolated temps
plt.title("Interpolated Temperatures 17070100")
plt.show()

代码运行没有错误，但我最终得到一张台湾的空地图。

我非常绝望，任何帮助将不胜感激！！！

4

1 回答 1

1

每当您向地图添加数据时，务必记住定义坐标系。在这种情况下，由于您的数据是纬度/经度，我将首先执行以下操作：

view.pcolormesh(..., transform=ccrs.PlateCarree())

您可能也有兴趣看到在使用 cartopy 绘制其他投影中的投影数据时认真使用的 transform 关键字。

于 2017-08-08T08:32:10.643 回答