Python - интерполяция температуры
В настоящее время я пытаюсь использовать код, предоставленный на этом веб-сайте ( https://unidata.github.io/MetPy/latest/examples/gridding/Point_Interpolation.html) для создайте карту Тайваня с линейной интерполяцией данных на ноутбуке 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()
Код выполняется без ошибок, но я получаю пустую карту Тайваня.
Я супер отчаянно, любая помощь будет высоко ценится!
1 ответ
Всякий раз, когда вы добавляете данные на картографическую карту, важно помнить, чтобы определить систему координат. В этом случае, так как ваши данные в широте / долготе, я бы начал делать что-то вроде:
view.pcolormesh(..., transform=ccrs.PlateCarree())
Вам также может быть интересно увидеть, как ключевое слово transform использовалось всерьез на графике спроецированных данных в других проектах с использованием картопий.