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

Maziy

王老师，我在您写的绘制雷达反射率的基础上进行修改，想计算出组合反射率，逻辑是对每一个仰角的rf进行垂直上的对比从而找出最大的反射率，但是在“km to degree“”这一块也采取同样的处理方法好像不太行，出来的图会花掉。请问该怎么处理呢？
测试脚本如下：
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 = 110.24989
lat =  19.99686

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

#scan从0-8计算出来的rf1-rf8分别是不同仰角下的基本反射率
#scan0
#Read data
fn = 'D:/180814-16radr/2018081413.32A'
f = addfile(fn)
scan = 0
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 j in xrange (len(e)):
    x[j,:] = dis * cos(e[j]) * cos(a[j])
    y[j,:] = dis * cos(e[j]) * sin(a[j])
x[nr,:] = x[0,:]
y[nr,:] = y[0,:]
rf1 = zeros((nr + 1, nd))
rf1[:nr,:nd] = rf
rf1[nr,:] = rf[0,:]

#scan1
#Read data
scan = 1
rf2  = f['reflectivity'][scan,:,:]
azi2 = f['azimuthR'][scan,:]
ele2 = f['elevationR'][scan,:]
#Get x/y (kilometers) coordinates of data
e2 = radians(ele2)
azi2 = 90 - azi2
a2 = radians(azi2)
nr2 = rf2.shape[0]
nd2 = rf2.shape[1]
x2 = zeros((nr2 + 1, nd2))
y2 = zeros((nr2 + 1, nd2))
for j in xrange (len(e2)):
    x2[j,:] = dis * cos(e2[j]) * cos(a2[j])
    y2[j,:] = dis * cos(e2[j]) * sin(a2[j])
x2[nr2,:] = x2[0,:]
y2[nr2,:] = y2[0,:]
rf12 = zeros((nr2 + 1, nd2))
rf12[:nr2,:nd2] = rf2
rf12[nr2,:] = rf2[0,:]
       
#将8个通道的基本反射率在垂直方向对比取大值，插值到地图上完成组合反射率       
test=rf1.copy()
xx=x.copy()
yy=y.copy()
#比较0和1
for i in range(0,361):
    for j in range(0,459):
        if        rf1[i,j]>rf12[i,j]:
            test[i,j]=rf1[i,j]
            xx[i,j]=x[i,j]
            yy[i,j]=y[i,j]
        elif rf1[i,j]==rf12[i,j]:
            test[i,j]=rf1[i,j]
            xx[i,j]=x[i,j]
            yy[i,j]=y[i,j]
        else:
            test[i,j]=rf12[i,j]
            xx[i,j]=x2[i,j]
            yy[i,j]=y2[i,j]
               


#km to degree
xx = xx * deg_1km + lon
yy = yy * deg_1km + lat

#Plot   
ax = axesm(bgcolor=[86,146,244])
#Add map layers
geoshow(19.37,110.10,size=8,color='r',marker='g')
text(109.95,19.42,u'屯昌',fontname=u'黑体',fontsize=20)
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=20, 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),(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)]
layer=pcolorm(xx, yy, test, 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(108,112) #控制绘图范围
ylim(17.7,20.7)
#xlim(lon - rd, lon + rd)
#ylim(lat - rd, lat + rd)
xaxis(tickin=False)
xaxis(tickvisible=False, location='top')
yaxis(tickin=False)
yaxis(tickvisible=False, location='right')
yticks(arange(18, 20.3,0.5))#控制出图y轴范围
#xticks(arange(108.3,111.3,1))
title('Radar reflectivity')

myj0219

感谢分享！！！！！！！！1

channam

请问这个图片是画的某个仰角的平面图还是？ 从第二段程序来看，x、y坐标是二维坐标，第一维代表了仰角，插值的时候是怎么操作呢？这一步没太看明白，rf,x,y = griddata((x,y), rf, method='nearest')，等号左边rf，x，y是什么意思？是对这三个要素都插值吗？从网上查到这个函数griddata(points, values, xi, method='linear', fill_value=numpy.nan, rescale=False)   xi是需要插值的空间，一般用 numpy.mgrid 函数生成后传入；我们这里不需要设置xi需要插值的空间吗？

channam

Maziy 发表于 2021-7-23 12:18
王老师，我在您写的绘制雷达反射率的基础上进行修改，想计算出组合反射率，逻辑是对每一个仰角的rf进行垂直 ...

请问这个图片是画的某个仰角的平面图还是？ 按理说（x，y）坐标应该有好几个不同仰角的值， 总看到把雷达数据画成这样经纬度的平面图，但仰角是怎么处理的呢？

Maziy

channam 发表于 2022-1-12 17:15
请问这个图片是画的某个仰角的平面图还是？ 按理说（x，y）坐标应该有好几个不同仰角的值， 总看到把雷达 ...

在上个版本的基础修改了，用了笨办法做循环画组合反射率

gaoyc

厉害