Python解析SWAN气象雷达数据--（解析、生成ASCII、Image、netCDF）
解析
from datetime import * import time import calendar import json import numpy as np from struct import * import binascii import netCDF4 file = open(r"D:/radarDataTest/Z_QPF_20140831203600.F030.bin", "rb") data = file.read(); print(len(data)) file.close() # file = open(r"D:/radarDataTest/Z_QPF_20140831203600.F030.bin", "rb") length = 0 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') length = length + 12+38+8+8 # print(zonName) print("数据说明: " + dataName) print("文件标志: " + flag) print("数据版本号: " + version) # year,month,day,hour,minute,interval, = unpack("HHHHHH", file.read(2+2+2+2+2+2)) print("时间: "+str(year)+"-"+str(month)+"-"+str(day)+" "+str(hour)+":"+str(minute)) print("时段长: "+str(interval)) length = length + 2+2+2+2+2+2 # XNumGrids,YNumGrids,ZNumGrids, = unpack("HHH", file.read(2+2+2)) print("X: " + str(XNumGrids)+" Y: "+str(YNumGrids)+" Z:"+str(ZNumGrids)) length = length + 2+2+2 # RadarCount, = unpack("i", file.read(4)) print("拼图雷达数: " + str(RadarCount)) length = length + 4 # 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)) length = length + 4+4+4+4+4+4 ZhighGrids=[] for i in range(0, 40): ZhighGrid, = unpack("f", file.read(4)) ZhighGrids.append(ZhighGrid) print("垂直方向的高度:"+str(ZhighGrids)) length = length + 40*4 # RadarStationNames=[] for i in range(0, 20): RadarStationName, = unpack("16s", file.read(16)) RadarStationName = RadarStationName.decode("gbk") RadarStationNames.append(RadarStationName.rstrip('\x00')) print("相关站点名称:"+str(RadarStationNames)) length = length + 20*16 # RadarLongitudes=[] for i in range(0, 20): RadarLongitude, = unpack("f", file.read(4)) RadarLongitudes.append(RadarLongitude) print("相关站点所在经度:"+str(RadarLongitudes)) length = length + 20*4 # RadarLatitudes=[] for i in range(0, 20): RadarLatitude, = unpack("f", file.read(4)) RadarLatitudes.append(RadarLatitude) print("相关站点所在纬度:"+str(RadarLatitudes)) length = length + 20*4 # RadarAltitudes=[] for i in range(0, 20): RadarAltitude, = unpack("f", file.read(4)) RadarAltitudes.append(RadarAltitude) print("相关站点所在海拔高度:"+str(RadarAltitudes)) length = length + 20*4 # MosaicFlags=[] for i in range(0, 20): MosaicFlag, = unpack("B", file.read(1)) MosaicFlags.append(MosaicFlag) print("该相关站点数据是否包含在本次拼图中:"+str(MosaicFlags)) length = length + 20*1 # m_iDataType, = unpack("h", file.read(2)) print("数据类型定义:"+str(m_iDataType)) if m_iDataType==0: print("unsigned char") elif m_iDataType==1: print("char") elif m_iDataType==2: print("unsigned short") elif m_iDataType==3: print("short") elif m_iDataType==4: print("unsigned short") length = length + 2 # m_iLevelDimension, = unpack("h", file.read(2)) print("每一层的向量数:"+str(m_iLevelDimension)) length = length + 2 # Reserveds=[] Reserveds, = unpack("168s", file.read(168)) Reserveds = Reserveds.decode("gbk").rstrip('\x00') print("该相关站点数据是否包含在本次拼图中: "+Reserveds) length = length + 168 #打印数据 valueZYX = [] for i in range(0, ZNumGrids): valueYX = [] for j in range(0, YNumGrids): valueX = [] for k in range(0, XNumGrids): value, = unpack("h", file.read(2)) #value, = unpack("b", file.read(1)) ''' if value > 0: print(value) ''' valueX.append(value) valueYX.append(valueX) valueZYX.append(valueYX) # #print("数据:"+str(valueZYX)) length = length + ZNumGrids*YNumGrids*XNumGrids*2 print(length) # print("----------------------------数据----------------------------") file.close()
生成ASCII
import time from struct import * 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) textZYX = [] for i in range(0, ZNumGrids): textYX = [] for j in range(0, YNumGrids): textX = [] for k in range(0, XNumGrids): value, = unpack("h", file.read(2)) textX.append(str(value)) textYX.append(' '.join(textX)) textZYX.append('\n'.join(textYX)) file.close() # #------------------------------------------------------------------------------- file_object = open('ASCIIData.txt', 'w') file_object.write("NCOLS " + str(XNumGrids) + "\n") file_object.write("NROWS " + str(YNumGrids) + "\n") file_object.write("XLLCENTER " + str(StartLon) + "\n") file_object.write("YLLCENTER " + str(StartLat - YReso * (YNumGrids - 1)) + "\n") # round(YReso, 3) * file_object.write("CELLSIZE " + str(XReso) + "\n") file_object.write("NODATA_VALUE " + str(-9999) + "\n") # # file_object.writelines(textZYX[0]) file_object.close() end = time.clock() print("read: %f s" % dateSpanTransfer) dateSpanTransfer = end - start #-------------------------------------------------------------------------------
生成Image（.img）
import time from struct import * from osgeo import gdal, osr from osgeo.gdalconst import * import numpy 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 in range(0, ZNumGrids): valueYX = [] for j in range(0, YNumGrids): valueX = [] for k in range(0, XNumGrids): value, = unpack("h", file.read(2)) valueX.append(value) valueYX.append(valueX) valueZYX.append(valueYX) file.close() # # #------------------------------------------------------------------------------- end = time.clock() dateSpanTransfer = end - start print("read: %f s" % dateSpanTransfer) # # driver = gdal.GetDriverByName('HFA') driver.Register() dataSetImg = driver.Create( "D:/radarDataTest/edarsImage.img", XNumGrids, YNumGrids, 1, gdal.GDT_Float32 ) # dataSetImg.SetGeoTransform( [ StartLon, XReso, 0, StartLat, 0, -YReso ] ) # srs = osr.SpatialReference() srs.SetWellKnownGeogCS( 'WGS84' ) dataSetImg.SetProjection( srs.ExportToWkt() ) # value2D = numpy.matrix( valueYX, dtype=numpy.float32 ) dataSetImg.GetRasterBand(1).WriteArray( value2D ) # dataSetImg = None #datasource.Destroy() #-------------------------------------------------------------------------------
生成netCDF
from datetime import * import time import calendar import json import numpy as np from struct import * import binascii import numpy from numpy.random import uniform from netCDF4 import Dataset 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 in range(0, ZNumGrids): valueYX = [] for j in range(0, YNumGrids): valueX = [] for k in range(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 in range(0, XNumGrids): for j in range(0, YNumGrids): for i in range(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 in range(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 in range(0, YNumGrids): y.append(StartLat - i * round(YReso, 3)) #y = numpy.array(y) lat[:] = y # z = [] for i in range(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()
