MeteoInfoLab脚本示例：绘制雷达反射率图

-rw-r----- 1 cliu cliu 2.6K Mar 30  2018 meteuse-grads.py
-rw-r----- 1 cliu cliu 4.0K Apr 10  2018 meteuse.py
-rw-r----- 1 cliu cliu  361 Mar 30  2018 piliang.py
-rw-r----- 1 cliu cliu  657 Mar 30  2018 plfile.py
-rw-r----- 1 cliu cliu  660 Mar 30  2018 pl.py
[cliu@c03n02 2017610nanjing]$ vi mete.py
[cliu@c03n02 2017610nanjing]$ vi meteuse.py
#import os.path
#import datetime
import os
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
def file_name(file_dir):
    for root, dirs, files in os.walk(file_dir):
        return files
datadir = 'data/'
fns = []
fns=file_name(datadir)

#Location of radar (Longitude/Latitude)
lon = 118.46
lat = 32.03

#Get degree per kilometer at the location
deg_1km = degree_1km(lon, lat)

#Read data
#fn = 'Z_RADR_I_Z9250_20170610084000_O_DOR_SA_CAP.bin'
for i in xrange (len(fns)):
#for i in xrange (8):
    fn=os.path.join(datadir,fns[i])
    print(fn)
    f=addfile(fn)
    rf=f['Reflectivity'][1,:,:]
    print(rf)
    azi=f['azimuthR'][0,:]
    dis=f['distanceR'][:]/1000.0
    ele=f['elevationR'][0,:]

    #Get x/y (kilometers) coordinates of data
    e=radians(ele)
    azi = 45 - azi
    a=radians(azi)
    x=zeros(rf.shape)
    y=zeros(rf.shape)
    for i in xrange (len(e)):
        x[i,:] = dis*cos(e[i])*cos(a[i])
        y[i,:] = dis*cos(e[i])*sin(a[i])

    #km to degree
    x = x * deg_1km + lon
    y = y * deg_1km + lat

    #Plot
    ##########clear memery####
    cla()
    ax = axesm(bgcolor=(204,255,255))
    #Add map layers
    #Set data folders
    mapdir = '/data/c03n02/cliu/software/meteinfo/MeteoInfo-1.4.8R7/MeteoInfo/map/CHN_adm'
    #Read shape files
    city_layer = shaperead(os.path.join(mapdir, 'CHN_adm2.shp'))
    #Plot
    geoshow(city_layer, edgecolor='lightgray',facecolor=(245,245,245),size=3)
    #geoshow(city_layer, edgecolor='lightgray',facecolor=(0,0,255),,size=6)
    city = geoshow('cn_cities', facecolor='r', size=28)
    city.addlabels('NAME', fontname=u'黑体', fontsize=36, yoffset=15)
    lab=city.getlabel(u'杭州')
    lab.getLegend().setText('Hangzhou')
    lab=city.getlabel(u'南京')
    lab.getLegend().setText('Nanjing')
    lab=city.getlabel(u'合肥')
    lab.getLegend().setText('Hefei')
    lab=city.getlabel(u'上海')
    lab.getLegend().setText('Shanghai')

    #geoshow('cn_province', edgecolor=None, facecolor=[230,230,230])
    #geoshow('cn_province', edgecolor=[80,80,80])
    #Add point
    #geoshow(15, -100, size=14, color='r', marker='S')
    #Add line
    #geoshow(lat, lon, size=2, color='b')
    #Add line and polygon
    #geoshow(lat, lon, color='r', size=2, linestyle=':')
    #geoshow(lat, lon, displaytype='polygon', color=[150,230,230,230], edgecolor='b', size=2)
    #Plot radar reflectivity
    levs=[5,10,15,20,25,30,35,40,45,50,55,60,65,70]
    #cols=[(255,255,255,0),(102,255,255),(0,162,232),(86,225,250),(3,207,14),\
    #    (26,152,7),(255,242,0),(217,172,113),(255,147,74),(255,0,0),\
    #    (204,0,0),(155,0,0),(236,21,236),(130,11,130),(184,108,208)]
    cols=[(255,255,255,0),(0,236,236),(1,160,246),(0,0,246),(0,255,0),\
        (0,200,0),(0,144,0),(255,255,0),(231,192,0),(255,144,0),\
        (255,0,0),(214,0,0),(192,0,0),(255,0,255),(153,85,201)]
    layer=pcolorm(x, y, rf, levs, colors=cols, order=1)
    colorbar(layer, shrink=0.8, label='dBZ', labelloc='top',fontsize=46)
    #Plot circles
    rr = array([50, 100, 150, 200])
    for r in rr:
        rd = r * deg_1km
        ax.add_circle((lon, lat), rd, edgecolor=(255,69,0),size=2)
    geoshow([lat,lat], [lon-rd,lon+rd], color=(255,69,0),size=3)
    geoshow([lat+rd,lat-rd], [lon,lon], color=(255,69,0),size=3)
    #Set plot
    #xlim(lon - rd, lon + rd)
    #ylim(lat - rd, lat + rd)
    xlim(117.968, 119.688)
    ylim(31.1842, 32.8038)
    xaxis(tickin=False)
    xaxis(tickvisible=False, location='top')
    yaxis(tickin=False)
    yaxis(tickvisible=False, location='right')
    xaxis(tickfontsize=46)
    yaxis(tickfontsize=46)
    yticks(arange(29, 34, 1))
    title('Radar reflectivity at '+fn[20:24]+'/'+fn[24:26]+'/'+fn[26:28]+'_'+fn[28:30]+':'+fn[30:32]+':'+fn[32:34],fontsize=46, fontname=u'楷体')
    #savefig('Z_RADR_2017-06-10_06_00.eps')
    savefig(fn[20:32]+'.png',3000,2000)