現在、このWebサイト(https://unidata.github.io/MetPy/latest/examples/gridding/Point_Interpolation.html#sphx-glr-examples-gridding-point-interpolation-py)で提供されているコードを使用して、Jupyterノートブックでデータの線形補間を行っている台湾の地図を作成しようとしています。Python - 温度補間
私のデータは、この形式である:
は#1707lat 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
私のコードは次のようになります:このフォームで
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
と、コードがエラーなしで実行
# 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()
しかし、私は台湾の空の地図に終わる。
私は非常に絶望的です、どんな助けも大歓迎です!
私はこのライブラリを知らないが、マップ自体の外にカラーマップを表現しているようだ。 'view.set_extent([120.5、122.5、24.5、25.5]) 'というコードは、座標を設定しているように見えます。しかし、マップの 'to_proj = ccrs.AlbersEqualArea(central_longitude = 120.0000、 central_latitude = 25.0000)の中心は、そのボックスの外側にあります... – agastalver