登录后查看更多精彩内容~
您需要 登录 才可以下载或查看，没有帐号？立即注册
x
本帖最后由 Masterpiece 于 2020-5-10 22:44 编辑
全球地形
程序中包含有matplotlib+basemap绘制contourf的基础操作，以及截取（truncate）部分官方色表（colormap）的做法，同时还有当contourf中设置extend='both'之后，重新定义低出绘制区间数据所使用的颜色的做法(set_under)。
- import numpy as np
- import matplotlib.pyplot as plt
- from mpl_toolkits.basemap import Basemap
- import netCDF4 as nc
- import matplotlib as mpl
- def truncate_colormap(cmap, minval=0.0, maxval=1.0, n=128):
- new_cmap = mpl.colors.LinearSegmentedColormap.from_list('trunc({n},{a:.2f},{b:.2f})'.format(n=cmap.name,a=minval, b=maxval),cmap(np.linspace(minval, maxval, n)))
- return new_cmap
- topo_file='ETOPO2v2g_f4.nc'
- data=nc.Dataset(topo_file)
- topo=data.variables['z'][:,:]
- lat=data.variables['y'][:]
- lon=data.variables['x'][:]
- lon_leftup=-180;lat_leftup=90
- lon_rightdown=180;lat_rightdown=-90 #建立地图投影
- fig, ax = plt.subplots(figsize=(14,9)) #建立绘图平台
- m = Basemap(projection='cyl', llcrnrlat=lat_rightdown, urcrnrlat=lat_leftup, llcrnrlon=lon_leftup, urcrnrlon=lon_rightdown, resolution='i')
- m.drawcoastlines(linewidth=0.2, color='gray',zorder=3)
- parallels = np.arange(-90,90+30,25) #纬线
- m.drawparallels(parallels,labels=[True,False,False,False],linewidth=0.01,dashes=[1,400])
- plt.yticks(parallels,len(parallels)*[''])
- meridians = np.arange(-180,180+30,30) #经线
- m.drawmeridians(meridians,labels=[False,False,False,True],linewidth=0.01,dashes=[1,400])
- plt.xticks(meridians,len(meridians)*[''])
- lons, lats = np.meshgrid(lon,lat) #经纬度2维化
- x, y = m(lons, lats) #投影映射
- cmap_new = truncate_colormap(plt.cm.terrain, 0.23, 1.0) #截取colormap，要绿色以上的（>=0.23）
- cmap_new.set_under([198/255,234/255,250/255]) #低于0的填色为海蓝
- lev=np.arange(0,5600,200)
- norm3 = mpl.colors.BoundaryNorm(lev, cmap_new.N) #标准化level，映射色标
- cf=m.contourf(x,y,topo,levels=lev,norm=norm3
- ,cmap=cmap_new,extend='both')
- cb=plt.colorbar(cf, ax=ax,shrink=0.7,aspect=30,pad=0.05,orientation='horizontal') #色标
- cb.ax.tick_params(labelsize=10,pad=2,direction='in') #色标tick
- plt.savefig('global_etopo.png',dpi=300,bbox_inches='tight') #存图
中国地形+气候160站
shp文件可以在家园里边搜（bou1_4l和country1），160站的csv在下边附件。
这个程序除了上边提到的技巧之外，其中一个就是地图投影的设置，还有使用汉字的操作。
还有就是使用shape文件maskout完美白化区域的做法，这个具体请参考下边给出的原贴。
matplotlib的scatter默认设置下的多边形marker，放大细看边角是圆的，假如利用path_effects为散点添加细边框可以解决这一问题，而且会让图更好看~
南海小图是通过在主图右下角添加图轴(axes)来达成的：plt.axes([小图x位置, 小图y位置, 小图的长, 小图的高])
- import numpy as np
- import matplotlib.pyplot as plt
- from mpl_toolkits.basemap import Basemap
- import netCDF4 as nc
- import matplotlib as mpl
- import pandas as pd
- import maskout
- import matplotlib.patheffects as path_effects
- def truncate_colormap(cmap, minval=0.0, maxval=1.0, n=128):
- new_cmap = mpl.colors.LinearSegmentedColormap.from_list('trunc({n},{a:.2f},{b:.2f})'.format(n=cmap.name,a=minval, b=maxval),cmap(np.linspace(minval, maxval, n)))
- return new_cmap
- topo_file='ETOPO2v2g_f4.nc'
- data=nc.Dataset(topo_file)
- print(data)
- topo=data.variables['z'][2600:,5400:]
- lat=data.variables['y'][2600:]
- lon=data.variables['x'][5400:]
- sta = pd.read_csv('sta160.csv', header=0, names=['sname', 'sid', 'slat', 'slon'])
- sta_name_160 = sta['sname']
- sta_id_160 = sta['sid']
- lo_160= sta['slon'].to_numpy()
- la_160 = sta['slat'].to_numpy()
- fig, ax = plt.subplots(figsize=(11,8)) #建立绘图平台
- m = Basemap(width=6500000,height=4300000
- ,resolution='l',projection='lcc',
- lat_1=30.,lat_2=45,lat_0=36,lon_0=109.)
- m.readshapefile('D:\\dt\\bou1_4l', 'chn', color='k', linewidth=0.5)
- m.drawcoastlines(linewidth=0.3, color='gray')
- parallels = np.arange(0,80,10) #纬线
- m.drawparallels(parallels,labels=[True,False,False,False],linewidth=0.3,dashes=[1,4])
- meridians = np.arange(80,150,10) #经线
- m.drawmeridians(meridians,labels=[False,False,False,True],linewidth=0.3,dashes=[1,4])
- lons, lats = np.meshgrid(lon,lat) #经纬度2维化
- x, y = m(lons, lats) #投影映射
- cmap_new = truncate_colormap(plt.cm.terrain, 0.23, 1.0) #截取colormap，要绿色以上的（>=0.23）
- cmap_new.set_under([198/255,234/255,250/255]) #低于0的填色为海蓝
- lev=np.arange(0,5600,200)
- norm3 = mpl.colors.BoundaryNorm(lev, cmap_new.N) #标准化level，映射色标
- cf=m.contourf(x,y,topo,levels=lev,norm=norm3
- ,cmap=cmap_new,extend='both')
- clip=maskout.shp2clip(cf,ax,m,'D:/dt/country1',['China'])
- cb=plt.colorbar(cf, ax=ax,shrink=0.7,aspect=30,pad=0.05,orientation='horizontal') #色标
- cb.ax.tick_params(labelsize=10,pad=2,direction='in') #色标tick
- aa,bb=m(lo_160, la_160)
- t=m.scatter(aa, bb, s=25, marker='^',label='CMA (%d Stations)'%len(lo_160), color='g', zorder=6)
- t.set_path_effects([path_effects.PathPatchEffect(edgecolor='k', facecolor='r', linewidth=0.3)])
- ax.legend(loc=1)
- # 南海小图
- a = plt.axes([0.73, 0.27, 0.12, 0.23])
- lon_leftup=107;lat_leftup=24
- lon_rightdown=121.3;lat_rightdown=2.4
- m = Basemap(projection='cyl', llcrnrlat=lat_rightdown, urcrnrlat=lat_leftup, llcrnrlon=lon_leftup, urcrnrlon=lon_rightdown, resolution='l')
- m.drawcoastlines(linewidth=0.3, color='gray')
- m.readshapefile('D:\\dt\\bou1_4l', 'chn', color='k', linewidth=0.5)
- x, y = m(lons, lats) #投影映射
- cf=m.contourf(x,y,topo,levels=lev,norm=norm3
- ,cmap=cmap_new,extend='both')
- clip2=maskout.shp2clip(cf,a,m,'D:/dt/country1',['China'])
- aa,bb=m(lo_160, la_160)
- t=m.scatter(aa, bb, s=25, marker='^',label='CMA (%d Stations)'%len(lo_160), color='g', zorder=6)
- t.set_path_effects([path_effects.PathPatchEffect(edgecolor='k', facecolor='r', linewidth=0.3)])
- font1={'family':'SimHei','size':8,'color':'black'}
- a.text(0.58,0.04,'南海诸岛',transform=a.transAxes,fontdict=font1,
- bbox=dict(boxstyle='square',ec='k',fc='w',pad=0.3))
- plt.savefig('chn_etopo160.png',dpi=300,bbox_inches='tight') #存图
- plt.show()
---中国气候160站地理位置信息---
sta160.csv
(4.23 KB, 下载次数: 116)
|
[复制链接]
窥视卡
雷达卡
发表于 2020-4-22 16:10:29
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
抢沙发
千斤顶
显身卡
发表于 2020-4-22 18:08:14
楼主