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

MeteoInfo · 发表于 2018-3-30 10:14:48

登录后查看更多精彩内容~

您需要登录才可以下载或查看，没有帐号？立即注册

x

本帖最后由 MeteoInfo 于 2019-5-16 14:22 编辑

CINRAD SA雷达数据读取和绘图，程序里有具体标注。

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 = 120.9586
lat = 31.0747
#Get degree per kilometer at the location
deg_1km = degree_1km(lon, lat)
#Read data
fn = 'D:/Temp/binary/Z_RADR_I_Z9002_20180304192700_O_DOR_SA_CAP.bin'
f = addfile(fn)
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 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,:]
#km to degree
x = x * deg_1km + lon
y = y * deg_1km + lat
#Plot
ax = axesm(bgcolor='b')
#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),(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(x, y, rf1, 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)
xaxis(tickin=False)
xaxis(tickvisible=False, location='top')
yaxis(tickin=False)
yaxis(tickvisible=False, location='right')
yticks(arange(29, 34, 1))
title('Radar reflectivity')

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 = 120.9586 
lat = 31.0747 
 
#Get degree per kilometer at the location 
deg_1km = degree_1km(lon, lat) 
 
#Read data 
fn = 'D:/Temp/binary/Z_RADR_I_Z9002_20180304192700_O_DOR_SA_CAP.bin' 
f = addfile(fn) 
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 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,:] 
 
#km to degree 
x = x * deg_1km + lon 
y = y * deg_1km + lat 
 
#Plot         
ax = axesm(bgcolor='b') 
#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),(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(x, y, rf1, 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) 
xaxis(tickin=False) 
xaxis(tickvisible=False, location='top') 
yaxis(tickin=False) 
yaxis(tickvisible=False, location='right') 
yticks(arange(29, 34, 1)) 
title('Radar reflectivity')

上述脚本中x, y坐标都是2维数组，形成的是非矩形网格，用pcolor/pcolorm函数形成的是多边形图集，也就是说每个非矩形网格是一个多边形（Polygon），这样做虽然无需进行数据插值，但网格太多的话会影响绘图速度。另一种数据处理方式是先将非矩形网格点（当作散点对待）插值为矩形网格，然后再用imshow/imshowm函数生成图像进行显示。插值会需要一些时间，但图像的显示会快很多。数据插值处理的脚本：

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 = 120.9586
lat = 31.0747
#Get degree per kilometer at the location
deg_1km = degree_1km(lon, lat)
#Read data
fn = 'D:/Temp/binary/Z_RADR_I_Z9002_20180304192700_O_DOR_SA_CAP.bin'
f=addfile(fn)
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(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
#Interpolate to griddata
rf,x,y = griddata((x,y), rf, method='nearest')
#rf,x,y = griddata((x, y), rf, method='idw', pointnum=20)
#rf,x,y = griddata((x, y), rf, method='idw', radius=0.1)
#Plot
ax = axesm(bgcolor='b')
#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),(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 = imshowm(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)
xaxis(tickin=False)
xaxis(tickvisible=False, location='top')
yaxis(tickin=False)
yaxis(tickvisible=False, location='right')
yticks(arange(29, 34, 1))
title('Radar reflectivity')

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 = 120.9586 
lat = 31.0747 
 
#Get degree per kilometer at the location 
deg_1km = degree_1km(lon, lat) 
 
#Read data 
fn = 'D:/Temp/binary/Z_RADR_I_Z9002_20180304192700_O_DOR_SA_CAP.bin' 
f=addfile(fn) 
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(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 
 
#Interpolate to griddata 
rf,x,y = griddata((x,y), rf, method='nearest') 
#rf,x,y = griddata((x, y), rf, method='idw', pointnum=20) 
#rf,x,y = griddata((x, y), rf, method='idw', radius=0.1) 
 
#Plot         
ax = axesm(bgcolor='b') 
#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),(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 = imshowm(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) 
xaxis(tickin=False) 
xaxis(tickvisible=False, location='top') 
yaxis(tickin=False) 
yaxis(tickvisible=False, location='right') 
yticks(arange(29, 34, 1)) 
title('Radar reflectivity')

♂雨已~○ · 发表于 2018-3-30 11:14:08

来了，顶！！！

finler · 发表于 2018-3-30 11:46:22

赞！！！

{:5_213:}{:5_213:}{:5_213:}
请问有径向速度没{:5_213:}{:5_213:}

chongzika · 发表于 2018-3-30 12:20:49

本帖最后由 chongzika 于 2018-3-30 12:25 编辑

王老师我用的1.4.8版本的报错，是版本问题么？
1
Running Jython script...
mipylib is loaded...
Traceback (most recent call last):
File "mete.py", line 30, in <module>
x[i,:] = dis*cos(e)*cos(a)
File "/data/c03n02/cliu/software/meteinfo/MeteoInfo-1.4.8/MeteoInfo/pylib/mipylib/numeric/miarray.py", line 221, in __setitem__
r = ArrayMath.setSection(self.array, ranges, value)
TypeError: setSection(): 3rd arg can't be coerced to ucar.ma2.Array, java.lang.Number

MeteoInfo · 发表于 2018-3-30 12:42:10

chongzika 发表于 2018-3-30 12:20
王老师我用的1.4.8版本的报错，是版本问题么？
1
Running Jython script...

需要用目前最新的版本

chongzika · 发表于 2018-3-30 15:05:12

MeteoInfo 发表于 2018-3-30 12:42
需要用目前最新的版本

#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
ax = axesm(bgcolor='b')
#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=15)
#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(x, y, rf, levs, colors=cols, order=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)
xaxis(tickin=False)
xaxis(tickvisible=False, location='top')
yaxis(tickin=False)
yaxis(tickvisible=False, location='right')
yticks(arange(29, 34, 1))
title('Radar reflectivity')
#savefig('Z_RADR_2017-06-10_06_00.eps')
savefig(fn[19:27]+'png',600,400)

复制代码

王老师，您好，我照着您说的对应雷达数据的脚本写了个循环，一下子处理200多个文件，但一般到处理10个左右就不行了内存。每次还需要释放内存么？我不太懂这个，具体错误请看
data/Z_RADR_I_Z9250_20170610003600_O_DOR_SA_CAP.bin
array([[-33.0, 35.5, 46.5, 45.0, 44.0, 40.5, 35.0, 34.5, 33.5, 35.0, 33.5, 31.5, 39.0, 30.5, 26.5, 35.0, 38.0, 34.0, 27.0, 19.5, 24.5, 30.5, 24.0, 23.5, 21.0, 20.0, 21.5, 21.5, 26.5, 26.5, 23.5, 30.0, 38.5, 46.0, 47.0, 45.0, 44.5, 46.5, 45.5, 44.0, 46.0, 39.0, 36.5, 35.5, 30.5, 31.5, 37.0, 43.0, 40.5, 45.0, 42.0, 34.5, 34.5, 34.5, 34.0, 34.5, 38.0, 34.5, 37.0, 38.0, 38.5, 41.5, 33.5, 27.5, 24.0, 28.5, 28.5, 31.5, 28.5, 30.5, 27.5, 29.5, 31.5, 30.5, 34.5, 37.0, 42.0, 43.5, 39.5, 41.5, 44.0, 38.5, 40.0, 40.0, 37.0, 40.0, 36.0, 39.5, 42.0, 44.0, 42.5, 41.0, 40.5, 40.5, 42.5, 41.5, 43.5, 45.0, 44.5, 45.5, 44.5, 46.5, 43.5, 39.0, 41.5, 41.0, 39.5, 41.5, 39.5, 38.5, 40.5, 39.0, 40.5, 37.0, 35.0, 38.5, 43.0, 39.0, 36.0, 31.0, 30.5, 35.0, 38.0, 38.5, 40.0, 39.5, 41.0, 43.0, 41.5, 40.5, 41.0, 44.0, 44.0, 44.5, 42.5, 45.5, 43.5, 42.0, 42.5, 40.0, 39.0, 42.0, 45.5, 45.5, 45.5, 45.0, 45.0, 48.0, 47.5, 45.5, 43.5, 39.5, 39.5, 40.5, 40.5, 40.0, 36.5, 34.0, 32.5, 34.0, 34.0, 33.5, 32.5, 34.0, 34.0, 34.0, 30.5, 32.0, 30.5, 29.5, 29.0, 28.0, 29.0, 27.0, 27.0, 26.5, 27.5, 27.0, 24.5, 24.5, 23.5, 23.0, 23.0, 22.5, 23.0, 26.0, 29.0, 30.0, 27.5, 26.5, 24.5, 27.0, 28.0, 31.0, 27.5, 27.5, 24.5, 27.0, 28.5, 26.0, 27.0, ...]])
Traceback (most recent call last):
 File "meteuse.py", line 63, in <module>
layer=pcolorm(x, y, rf, levs, colors=cols, order=1)
 File "/data/c03n02/cliu/software/meteinfo/MeteoInfo-1.4.8R7/MeteoInfo/pylib/mipylib/plotlib/miplot.py", line 5140, in pcolorm
return gca.pcolor(*args, **kwargs)
 File "/data/c03n02/cliu/software/meteinfo/MeteoInfo-1.4.8R7/MeteoInfo/pylib/mipylib/plotlib/axes.py", line 540, in pcolor
layer = ArrayUtil.meshLayer(x.asarray(), y.asarray(), a.asarray(), ls, lonlim)
 at org.meteoinfo.shape.PolygonShape.updatePolygons(PolygonShape.java:338)
 at org.meteoinfo.shape.PolygonShape.setPoints(PolygonShape.java:188)
 at org.meteoinfo.data.ArrayUtil.meshLayer(ArrayUtil.java:1907)
 at sun.reflect.NativeMethodAccessorImpl.invoke0(Native Method)
 at sun.reflect.NativeMethodAccessorImpl.invoke(NativeMethodAccessorImpl.java:57)
 at sun.reflect.DelegatingMethodAccessorImpl.invoke(DelegatingMethodAccessorImpl.java:43)
 at java.lang.reflect.Method.invoke(Method.java:606)
java.lang.OutOfMemoryError: java.lang.OutOfMemoryError: GC overhead limit exceeded

MeteoInfo · 发表于 2018-3-30 16:05:29

chongzika 发表于 2018-3-30 15:05
王老师，您好，我照着您说的对应雷达数据的脚本写了个循环，一下子处理200多个文件，但一般到处理10个左 ...

在51行后面加上cla()试试。cla是清除之前的axes
#Plot
cla()

chongzika · 发表于 2018-3-30 18:20:42

MeteoInfo 发表于 2018-3-30 16:05
在51行后面加上cla()试试。cla是清除之前的axes
#Plot
cla()

可以了，多谢王老师

chongzika · 发表于 2018-3-30 18:51:38

MeteoInfo 发表于 2018-3-30 16:05
在51行后面加上cla()试试。cla是清除之前的axes
#Plot
cla()

顺便问一下王老师，图中的汉字是自动绘制的，能不能改成汉语拼音呢

freekiller · 发表于 2018-3-31 10:51:17

厉害了

		自动登录	找回密码
密码			立即注册