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

独孤酒见

厉害了，王老师。

liuyg123

太棒了 感谢王老师

zyl1990

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 = 108.73
lat = 19.22

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

#Read data
fn = 'D:/dongfang_rad/Z_RADR_I_Z9072_20180722004200_O_DOR_SA_CAP.bin'
f=addfile(fn)
rf=f['Reflectivity'][1,:,:]
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 j in xrange (len(e)):
    x[j,:] = dis * cos(e[j]) * cos(a[j])
    y[j,:] = dis * cos(e[j]) * sin(a[j])

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

#Plot        
ax = axesm(bgcolor='w')
#Add map layers
#geoshow('cn_province', edgecolor=None, facecolor=[230,230,230])
geoshow('cn_province', edgecolor=[80,80,80])
#city = geoshow('cn_cities', facecolor='r', size=8)
#city.addlabels('NAME', fontname=u'黑体', fontsize=16, yoffset=18)
#Plot radar reflectivity
levs=[5,10,15,20,25,30,35,40,45,50,55,60,65,70]
cols=[(255,255,255,0),(0,255,255),(0,157,255),(0,0,255),(9,130,175),\
    (0,255,0),(8,175,20),(255,214,0),(255,152,0),(255,0,0),\
    (221,0,27),(188,0,54),(121,0,109),(121,51,160),(195,163,212)]
layer=pcolorm(x, y, rf, levs, colors=cols, zorder=1)
colorbar(layer, shrink=0.8, 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')
#Set plot
#xlim(lon - rd, lon + rd)
#ylim(lat - rd, lat + rd)
xlim(108, 111.5)
ylim(18, 20.5)
xaxis(tickin=False)
xaxis(tickvisible=False, location='top')
yaxis(tickin=False)
yaxis(tickvisible=False, location='right')
yticks(arange(18, 20.5, 1))
title('Radar reflectivity')
王老师，按照您的脚本画出来的反射率因子，为什么总会出现锥形空白呢？

MeteoInfo

zyl1990 发表于 2018-9-17 20:43
def degree_1km(lon, lat):
    '''
    Get degree per kilometer at a location.   

1楼原脚本程序有些问题，已经更新，你再试试。

zyl1990

MeteoInfo 发表于 2018-9-18 09:26
1楼原脚本程序有些问题，已经更新，你再试试。

谢谢王老师，有组合反射率脚本吗

zyl1990

MeteoInfo 发表于 2018-9-18 09:26
1楼原脚本程序有些问题，已经更新，你再试试。

谢谢王老师，有组合反射率脚本吗

zyl1990

MeteoInfo 发表于 2018-9-18 09:26
1楼原脚本程序有些问题，已经更新，你再试试。

谢谢王老师，有组合反射率脚本吗

MeteoInfo

zyl1990 发表于 2018-9-18 11:06
谢谢王老师，有组合反射率脚本吗

自己先尝试吧

MeteoInfo

zyl1990 发表于 2018-9-18 11:06
谢谢王老师，有组合反射率脚本吗

自己先尝试吧

MeteoInfo

zyl1990 发表于 2018-9-18 11:06
谢谢王老师，有组合反射率脚本吗

自己先尝试吧