2017-08-08 16 views
-1

現在、この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() 

しかし、私は台湾の空の地図に終わる。

私は非常に絶望的です、どんな助けも大歓迎です!

+0

私はこのライブラリを知らないが、マップ自体の外にカラーマップを表現しているようだ。 'view.set_extent([120.5、122.5、24.5、25.5]) 'というコードは、座標を設定しているように見えます。しかし、マップの 'to_proj = ccrs.AlbersEqualArea(central_longitude = 120.0000、 central_latitude = 25.0000)の中心は、そのボックスの外側にあります... – agastalver

答えて

0

カートリストマップにデータを追加するときは常に、座標系を定義することを忘れないようにしてください。あなたはまた、Plotting projected data in other projectons using cartopyで本格的に使用されている変換キーワードを見ることに興味があるかもしれません

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

:あなたのデータは、緯度/経度であるので、この場合、私のような何かをすることから始めたいです。