2014-08-27 21 views
0

我正在尝试使用streamplot函数来绘制带有底图的风场,投影"ortho"。我的测试代码主要是基于下面这个例子: Plotting wind vectors and wind barbsstreamplot不能与matplotlib底图一起使用

这里是我的代码:

import numpy as np 
import matplotlib.pyplot as plt 
import datetime 
from mpl_toolkits.basemap import Basemap, shiftgrid 
from Scientific.IO.NetCDF import NetCDFFile as Dataset 

# specify date to plot. 
yyyy=1993; mm=03; dd=14; hh=00 
date = datetime.datetime(yyyy,mm,dd,hh) 
# set OpenDAP server URL. 
URLbase="http://nomads.ncdc.noaa.gov/thredds/dodsC/modeldata/cmd_pgbh/" 
URL=URLbase+"%04i/%04i%02i/%04i%02i%02i/pgbh00.gdas.%04i%02i%02i%02i.grb2" %\ 
     (yyyy,yyyy,mm,yyyy,mm,dd,yyyy,mm,dd,hh) 
data = Dataset(URL) 
#data = netcdf.netcdf_file(URL) 
# read lats,lons 
# reverse latitudes so they go from south to north. 
latitudes = data.variables['lat'][:][::-1] 
longitudes = data.variables['lon'][:].tolist() 
# get wind data 
uin = data.variables['U-component_of_wind_height_above_ground'][:].squeeze() 
vin = data.variables['V-component_of_wind_height_above_ground'][:].squeeze() 
# add cyclic points manually (could use addcyclic function) 
u = np.zeros((uin.shape[0],uin.shape[1]+1),np.float64) 
u[:,0:-1] = uin[::-1]; u[:,-1] = uin[::-1,0] 
v = np.zeros((vin.shape[0],vin.shape[1]+1),np.float64) 
v[:,0:-1] = vin[::-1]; v[:,-1] = vin[::-1,0] 
longitudes.append(360.); longitudes = np.array(longitudes) 
# make 2-d grid of lons, lats 
lons, lats = np.meshgrid(longitudes,latitudes) 
# make orthographic basemap. 
m = Basemap(resolution='c',projection='ortho',lat_0=60.,lon_0=-60.) 
# create figure, add axes 
fig1 = plt.figure(figsize=(8,10)) 
ax = fig1.add_axes([0.1,0.1,0.8,0.8]) 
# define parallels and meridians to draw. 
parallels = np.arange(-80.,90,20.) 
meridians = np.arange(0.,360.,20.) 
# first, shift grid so it goes from -180 to 180 (instead of 0 to 360 
# in longitude). Otherwise, interpolation is messed up. 
ugrid,newlons = shiftgrid(180.,u,longitudes,start=False) 
vgrid,newlons = shiftgrid(180.,v,longitudes,start=False) 
# now plot. 
lonn, latt = np.meshgrid(newlons, latitudes) 
x, y = m(lonn, latt) 
st = plt.streamplot(x, y, ugrid, vgrid, color='r', latlon='True') 
# draw coastlines, parallels, meridians. 
m.drawcoastlines(linewidth=1.5) 
m.drawparallels(parallels) 
m.drawmeridians(meridians) 
# set plot title 
ax.set_title('SLP and Wind Vectors '+str(date)) 
plt.show() 

后运行的代码,我得到了一个空白的地图,在左下角的红色涂抹(请参阅该图)。将这些污迹放大后,我可以在平面投影中看到风流(而不是“正射影”)。所以我猜这是地图上数据投影的问题。我曾尝试功能transform_vector,但它不能解决问题enter image description here有谁可以告诉我,我做错了什么,请!谢谢。

更新代码之后的新地图: enter image description here

回答

1

您正在密谋lat/lon协调与正交投影在地图上。通常情况下,你可以通过改变你的绘图命令来解决这个问题:

m.streamplot(mapx, mapy, ugrid, vgrid, color='r', latlon=True) 

但是,你的坐标阵列不具有相同的尺寸,需要进行固定的。

+0

谢谢@Rutger Kassies。我添加了两行代码:'loon,latt = np.meshgrid(newlines,latitudes)x,y = m。(lonn,latt)',然后'm.streamplot(x,y,ugrid,vgrid,color =' r',latlon = True)'。现在我得到一个空白的地图,其中包含一个奇怪的行,并在屏幕上显示一些消息:'package/numpy/ma/core.py:802:RuntimeWarning:遇到的值无效 return umath.less(x,self.critical_value)''。请参阅上面的更新代码和新的空白地图 – 2014-08-27 15:48:07

相关问题