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()