- 积分
- 8083
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2011-12-27
- 最后登录
- 1970-1-1
|
5金钱
本帖最后由 小其其格 于 2019-6-10 08:26 编辑
- 文章有些图是用MeteoInfo绘制的,最后见刊排版阶段,编辑部要求X坐标的经度数值只在最后
- 显示°E,Y坐标只在最后数值显示°N。按照MeteoInfo官网对xticks函数的说明修改的:
xticks(arange(118.5, 122, 1.0), ['118.5','119.5','120.5','121.5E'],fontsize=fontSize, fontname='Times New Roman')
@MeteoInfo
绘图代码:
- dataDir = 'E:/wfl/MI/fig2/a~b/'
- fileType = '.bin'
- #字体大小
- fontSize = 36
- import os,sys,datetime
- meteoInfoLabName = os.path.abspath(sys.argv[0])
- print(meteoInfoLabName)
- #读取指定文件夹中的指定文件类型的文件名
- def get_filename(path,fileType):
- file_name =[]
- final_name = []
- for files in os.listdir(path): #root为目录路径 #dirs为路径下的子目录 #files为路径下的所有非子目录
- if fileType in files:
- file_name.append(files.replace(fileType,''))#生成不带‘.bin’后缀的文件名组成的列表
- final_name = [path +item +fileType for item in file_name]#生成‘.bin’后缀的文件名组成的绝对路径列表
- return final_name,file_name #输出列表
- def degree_1km(lon, lat):
- '''
- Get degree per kilometer at a location.
- '''
- dis_km = distance([lon,lon+1], [lat,lat], islonlat=True) / 1000
- return 1 / dis_km
- #Location of radar (Longitude/Latitude)
- lon = 121.514
- lat = 30.070
- minlat=28.0
- minlon=118.5
- maxlat=30.5
- maxlon=122.0
- #Get degree per kilometer at the location
- deg_1km = degree_1km(lon, lat)
- #读取多个文件的绝对路径和文件名
- fns,files = get_filename(dataDir,fileType) #fns为带.bin的文件名的文件列表
- nFiles = len(fns) #文件数量
- print(nFiles)
- print(fns)
- print('\n')
- print(files)
- #读取城市经纬度和名称等信息
- #Read city data of 11 stations from data file
- datafn = dataDir+'CityVars.Zhejiang.txt'
- ncol = numasciicol(datafn)
- nrow = numasciirow(datafn)
- var = asciiread(datafn, shape=(nrow,ncol))
- var_city = var[:,2]
- #Read station name and lon/lat
- stfn = dataDir+'StationNames.Zhejiang.csv'
- table = readtable(stfn, delimiter=',', format='%i%2s%2f',encoding='gb2312')
- stnames = table['Name']
- lat_city = table['LAT']
- lon_city = table['LON']
- #Read Radar datas and plot
- for k in range(nFiles):
-
- f = addfile(fns[k])
- scan = 1
- rf = f['Reflectivity'][scan,:,:]
- azi = f['azimuthR'][scan,:]
- dis = f['distanceR'][:]/1000.0
- ele = f['elevationR'][scan,:]
- #Get x/y (kilometers) coordinates of data
- e = radians(ele)
- azi = 90 - azi
- a = radians(azi)
- nr = rf.shape[0]
- nd = rf.shape[1]
- x = zeros((nr + 1, nd))
- y = zeros((nr + 1, nd))
- for i in xrange (len(e)):
- x[i,:] = dis * cos(e[i]) * cos(a[i])
- y[i,:] = dis * cos(e[i]) * sin(a[i])
- x[nr,:] = x[0,:]
- y[nr,:] = y[0,:]
- rf1 = zeros((nr + 1, nd))
- rf1[:nr,:nd] = rf
- rf1[nr,:] = rf[0,:]
- #km to degree
- x = x * deg_1km + lon
- y = y * deg_1km + lat
- #Plot
- #ax = axesm(bgcolor='b')
- ax = axesm()
- #Add map layers
- geoshow('cn_province', edgecolor=None, facecolor=[230,230,230])
- geoshow('cn_province', edgecolor=[80,80,80])
-
- #Plot radar reflectivity
- levs=[0,5,10,15,20,25,30,35,40,45,50,55,60,65]
- #color bar of CINRAD
- cols=[(255,255,255,0),(0,236,238),(0,161,244),(0,0,246),(0,255,0),\
- (0,200,0),(0,145,0),(255,255,0),(230,193,0),(254,146,0),\
- (254,0,0),(215,0,0),(192,0,0),(255,0,254),(155,85,199)]
- layer=pcolorm(x, y, rf1, levs, colors=cols, zorder=1)
- colorbar(layer, orientation='horizontal',shrink=0.8,fontsize=28,fontname=u'Times New Roman',label='dBZ', labelloc='top')
- #Plot circles
- rr = array([50, 100, 150, 200])
- for r in rr:
- rd = r * deg_1km
- ax.add_circle((lon, lat), rd, edgecolor='r')
- geoshow([lat,lat], [lon-rd,lon+rd], color='r')
- geoshow([lat+rd,lat-rd], [lon,lon], color='r')
- #绘制城市名称和描点
- layer_city = scatterm(lon_city, lat_city, var_city, size=5)
- layer_city.addfield('Name', 'string', stnames)
- layer_city.addlabels('Name', fontname=u'楷体', fontsize=fontSize, yoffset=-8)
-
- #Set plot
- if k ==3:
- minlat=28.0
- minlon=120.0
- maxlat=30.5
- maxlon=123.5
- xaxis(linewidth = 1.6) #坐标轴粗细
- yaxis(linewidth = 1.6) #坐标轴粗细
- #xlim(minlon, maxlon)
- ylim(minlat, maxlat)
- yticks(arange(minlat, maxlat,0.5), fontsize=fontSize, fontname='Times New Roman') #arange( start, stop, step )
- xticks(arange(118.5, 122, 1.0), ['118.5','119.5','120.5','121.5E'],fontsize=fontSize, fontname='Times New Roman')
-
- #保存图像
- #savefig(dataDir+files[k]+'.png', width=1920, height=1080, dpi=600)
- savefig(dataDir+files[k]+'.eps', width=1920, height=1080, dpi=600)
- print('Drawing ' +files[k] +'done!')
- if k < nFiles-1:
- cla() #清除,不然title要重叠
-
复制代码
出来的图如下,X坐标仍旧是每个数值后面有°E……
|
|