关于Basemap绘制天气图的问题

po_po1

请问大神们，现有一个全球的nc数据，例如：GTOPO30_10min.nc，该数据表征全球地形。
我只想绘制中国区域的地形，该进行哪些操作，才能只绘出中国区域地形？

po_po1

我只能绘出全球地形图：
程序：
from netCDF4 import Dataset
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
q=Dataset(r'C:/Users/Administrator/Desktop/GTOPO30_10min.nc')  #数据在桌面放着
a=q.variables['HT'][:]                                        #读取的地形数据
fig=plt.figure(figsize=(8,8))
m=Basemap(projection='cyl')
m.drawcountries()
m.drawcoastlines(linewidth=0.4)
ny=a.shape[0]
nx=a.shape[1]
x,y=m.makegrid(nx,ny)
lons,lats=m(x,y)
cs=plt.contourf(lons,lats,a)
plt.show()
q.close()
#######################################
#######################################
#######################################
#######################################
2014年7月24修改
最终代码：
from netCDF4 import Dataset
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
q=Dataset(r'D:/Desktop/GTOPO30_2min.nc')
a=q.variables['HT'][:]
lat=q.variables['lat'][:]
lon=q.variables['lon'][:]
fig=plt.figure(figsize=(8,8))
fig.add_axes([0.05,0.05,0.9,0.9])
m=Basemap(projection='cyl',llcrnrlon=70, urcrnrlon=135,\
          llcrnrlat=10, urcrnrlat=56, lon_0=125, lat_0=25, resolution='l')
m.drawcountries(color='w')
m.drawcoastlines(linewidth=0.4,color='w')
m.drawparallels(circles=np.arange(-80,90,20),labels=[0,1,0,0])
m.drawmeridians(meridians=np.arange(0,360,20),lables=[0,0,0,1],fontsize=10)
ny=a.shape[0]
nx=a.shape[1]
map0=r'D:\Program Files\WinPython-64bit-2.7.5.3\data\China boundary\province boundary\bou2_4l'
m.readshapefile(map0,'world',linewidth=0.5)
########
########设定绘图区域
index1=np.logical_and(lat>=10,lat<=56)#数组与常数逻辑运算（比较大小），得到布尔数组
index2=np.logical_and(lon>=70,lon<=135)
y=lat[index1]
x=lon[index2]
a=a[index1,:]
a=a[:,index2]
########
cs=plt.contourf(x,y,a)
cbar=m.colorbar(cs,location='bottom',pad=0.1,extend='both',extendfrac='auto')
plt.show()
q.close()

river

请问您这数据哪里下载的，我也去下载一个

阿阿飞飞

修改m=Basemap(projection='cyl')为：
map = Basemap(projection='cyl',llcrnrlon=100, urcrnrlon=150, llcrnrlat=0, urcrnrlat=50, lon_0=125, lat_0=25, resolution='l')
自己根据需要的经纬度范围进行调整，上面参数项含义为：起始和截止经度；起始和截止纬度；中心经度；中心纬度，分辨率

kongfeng0824

不太懂python，帮顶了

po_po1

阿阿飞飞 发表于 2014-2-19 12:41

                               
登录/注册后可看大图

修改m=Basemap(projection='cyl')为：
map = Basemap(projection='cyl',llcrnrlon=100, urcrnrlon=150, ll ...

不行啊，这样只是把地图画成中国区域了，数据还是全球的，成这个样子了：

po_po1

river 发表于 2014-2-19 12:32

                               
登录/注册后可看大图

请问您这数据哪里下载的，我也去下载一个

模式运行所需数据，网址：
http://users.ictp.it/~pubregcm/RegCM4/globedat.htm#part4
地形数据在第7组的GTOPO30 Topography下

阿阿飞飞

同时要对你读取出来的经纬度进行分割，在x,y=m.makegrid(nx,ny)语句后面添加：
index=np.logical_and(np.logical_and(x>100,x<150),np.logical_and(y>0,y<50))
x=x[index]
y=y[index]
a=a[index]

river

po_po1 发表于 2014-2-19 14:08

                               
登录/注册后可看大图

模式运行所需数据，网址：
http://users.ictp.it/~pubregcm/RegCM4/globedat.htm#part4
地形数据在第7组 ...

谢谢，我已经下载了。你可以试一下meteoinfo这个软件，可以利用它屏蔽中国外的数据

po_po1

river 发表于 2014-2-19 19:27

                               
登录/注册后可看大图

谢谢，我已经下载了。你可以试一下meteoinfo这个软件，可以利用它屏蔽中国外的数据

meteoinfo性能和GrADs类似吗？好不好学啊？