I am currently trying to use the code provided on this website (https://unidata.github.io/MetPy/latest/examples/gridding/Point_Interpolation.html#sphx-glr-examples-gridding-point-interpolation-py) to create a map of taiwan with the linear interpolation of data on a Jupyter notebook.
My data is in this form:
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
and also in this form:
#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
My code looks like this:
# 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()
The code runs without an error, but I end up with an empty map of Taiwan.
I am super desperate, any help would be much appreciated!!!
