python解析雷达数据_Python解析SWAN气象雷达数据--（解析、生成ASCII、Image、netCDF）...

from datetime import *

importtimeimportcalendarimportjsonimportnumpy as npfrom struct import *

importbinasciiimportnumpyfrom numpy.random importuniformfrom netCDF4 importDataset

start=time.clock()

file= open(r"D:/radarDataTest/Z_QPF_20140831203600.F030.bin", "rb")#zonName,dataName,flag,version, = unpack("12s38s8s8s", file.read(12+38+8+8))

zonName= zonName.decode("gbk").rstrip('\x00')

dataName= dataName.decode("gbk").rstrip('\x00')

flag= flag.decode("gbk").rstrip('\x00')

version= version.decode("gbk").rstrip('\x00')#print(zonName)print("数据说明:" +dataName)print("文件标志:" +flag)print("数据版本号:" +version)#length =0

length= length + 2+2+2+2+2+2 #时间说明

file.read(length)

XNumGrids,YNumGrids,ZNumGrids,= unpack("HHH", file.read(2+2+2))print("X:" + str(XNumGrids)+"Y:"+str(YNumGrids)+"Z:"+str(ZNumGrids))

length=0

length= length + 4 #拼图雷达数

file.read(length)#StartLon,StartLat,CenterLon,CenterLat,XReso,YReso, = unpack("ffffff", file.read(4+4+4+4+4+4))print("开始经度:"+str(StartLon)+"开始纬度:"+str(StartLat)+"中心经度:"+str(CenterLon)+"中心纬度:"+str(CenterLat))print("经度方向分辨率:"+str(XReso)+"纬度方向分辨率:"+str(YReso))

ZhighGrids=[]for i in range(0, 40):

ZhighGrid,= unpack("f", file.read(4))

ZhighGrids.append(ZhighGrid)print("垂直方向的高度:"+str(ZhighGrids))#length =0

length= length + 20*16 #相关站点名称

length = length + 20*4 #相关站点所在经度

length = length + 20*4 #相关站点所在纬度

length = length + 20*4 #相关站点所在海拔高度

length = length + 20*1 #该相关站点数据是否包含在本次拼图中

length = length + 2 #数据类型定义

length = length + 2 #每一层的向量数

length = length + 168 #保留信息

file.read(length)

valueZYX=[]for i inrange(0, ZNumGrids):

valueYX=[]for j inrange(0, YNumGrids):

valueX=[]for k inrange(0, XNumGrids):#value, = unpack("h", file.read(2))

#textX.append(str(value/10.0))

value, = unpack("b", file.read(1))

textX.append(str(value*2+66.0))

valueYX.append(valueX)

valueZYX.append(valueYX)

file.close()#valueXYZ =[]for k inrange(0, XNumGrids):for j inrange(0, YNumGrids):for i inrange(0, ZNumGrids):

valueXYZ.append(valueZYX[i][j][k])#file = open(r"D:/radarDataTest/Z_QPF_20140831203600.F030.bin", "rb")

rootgrp= Dataset("test.nc", "w", format="NETCDF4")#rootgrp = Dataset("test.nc", "a")#fcstgrp = rootgrp.createGroup("forecasts")

lon= rootgrp.createDimension("lon", XNumGrids)

lat= rootgrp.createDimension("lat", YNumGrids)

alt= rootgrp.createDimension("alt", ZNumGrids)

lon= rootgrp.createVariable("lon", "f8", ("lon",))

lat= rootgrp.createVariable("lat", "f8", ("lat",))

alt= rootgrp.createVariable("alt", "f8", ("alt",))#val = rootgrp.createVariable("val","f4",("zz","yy","xx",))

val = rootgrp.createVariable("val","f4",("lon","lat","alt",))#rootgrp.description =dataName

rootgrp.history= "创建时间：" + time.strftime('%Y-%m-%d %X', time.localtime())

rootgrp.Source_Software= "SmartMap"

#lon.units = "degrees_east"lon.long_name= "longitude coordinate"lon.standard_name= "longitude"

#lat.units = "degrees_north"lat.long_name= "latitude coordinate"lat.standard_name= "latitude"

#alt.units = "m"alt.long_name= "altitude"alt.standard_name= "heigh"

#val.long_name = "value"val.esri_pe_string= 'GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]'val.coordinates= "lon lat alt"val.units= "Degree"val.missing_value= -9999

#interval = 0.009999999776482582

interval = 0.01

#x = numpy.arange(-90,91,2.5)

x=[]for i inrange(0, XNumGrids):

x.append(StartLon+ i * round(XReso, 3))#x = numpy.array(x)

lon[:] =x#

#y = numpy.arange(-180,180,2.5)

y =[]for i inrange(0, YNumGrids):

y.append(StartLat- i * round(YReso, 3))#y = numpy.array(y)

lat[:] =y#z=[]for i inrange(0, ZNumGrids):

z.append(ZhighGrids[i])#z = numpy.array(z)

alt[:] =z#

#kk = uniform(size=(2,3,4,5))#print(kk)

#val[::]=valueZYX

val[::] =valueXYZ#rootgrp.close()