PERSIANN 降雨数据使用教程

一、前言

PERSIANN，“使用人工神经网络从遥感信息中估算降水”，是一种基于卫星的降水检索算法，可提供近乎实时的降雨信息。

该算法使用来自全球地球同步卫星的红外 (IR) 卫星数据作为降水信息的主要来源。 红外图像的降水基于云顶温度和降水率之间的统计关系。 然后使用低地球轨道卫星（例如，热带降雨测量任务微波成像仪、特殊传感器微波成像仪、高级微波扫描辐射计-地球观测系统）提供的卫星微波数据校准基于红外线的降水估计。 校准技术依赖于一种自适应训练算法，该算法在微波观测可用时（大约每隔 3 小时）更新检索参数。

PERSIANN 数据可以二进制格式获取。 这种格式不同于我们之前看到的geotif、netCDF和HDF。 因此，需要一种不同的方法来提取数据

二、数据下载

转到 ftp 站点并下载 2000 年 3 月的数据。

ftp://persiann.eng.uci.edu/pub/PERSIANN/tar_6hr/

提取并解压文件夹 d:/persiann/ 中的文件.

当使用 tar 提取然后解压缩时，6 小时文件将具有以下名称：raw6hrYYDOYHS.bin

YY : 2 位数年份，从 00 开始表示 2000

DOY：一年中的某一天

HS：6小时累计周期的起始时间，(0,6,12,18)

例如，文件 raw6hr0512518.bin 表示 2005 年 5 月 5 日 18 到 24 之间累积的降水量

数据以 C 样式以行为中心的格式存储。 第一个值以 49.875、0.125 为中心，第二个值以 49.875、0.375 为中心，最后两个值以：-49.875、359.625 和 -49.875、359.875 为中心

数据采用来自 SUN 系统（大端）的 4 字节二进制浮点数。从二进制文件中提取二进制数据需要执行以下步骤。 它基本上逐行扫描数据并将数据放入矩阵中。

三、数据处理

（1）在 d:\training1\scripts\ 中创建脚本 script5a.py 并运行下面的代码。

# import libraries
import matplotlib.pyplot as plt
import numpy as np
from struct import unpack
from osgeo import gdal
 
# inputfile
BinaryFile = "/path/to/raw6hr0006100.bin" # '3B42_daily.2009.05.31.7.bin'. Make sure you adjust the location.
 
# open binary file
f = open(BinaryFile, "rb")
 
# set file dimensions
xs = 1440
ys = 400
 
# set number of bytes in file
NumbytesFile = xs * ys
 
# number of columns in row
NumElementxRecord = -xs
 
# create empty array to put data in
myarr = []
 
# loop trough the binary file row by row
for PositionByte in range(NumbytesFile,0, NumElementxRecord):
 
        Record = ''
 
        # the dataset starts at 0 degrees, use 720 to convert to -180 degrees
        for c in range (PositionByte-720, PositionByte, 1):
                f.seek(c * 4)
                DataElement = unpack('>f', f.read(4))
                Record = Record  + str("%.2f" % DataElement + ' ')
 
        # 0 - 180 degrees
        for c in range (PositionByte-1440 , PositionByte-720, 1):
                f.seek(c * 4)
                DataElement = unpack('>f', f.read(4))
                Record = Record  + str("%.2f" % DataElement + ' ')
 
        # add data to array
        myarr.append(Record[:-1].split(" "))
 
# close binary file
f.close()
 
# Array to numpy float
myarr = np.array(myarr).astype('float')
 
# mirror array
myarr = myarr[::-1]
 
# show data
plt.imshow(myarr)
plt.clim(0,10)
plt.show()

如果出现带有降雨的图像，则您已成功将二进制数据导入到 numpy 数组中。

（2）在 d:\training1\scripts\ 中创建脚本 script5b.py 并运行下面的代码。 确保更正输入和输出路径。

# import libraries
import matplotlib.pyplot as plt
import numpy as np
from struct import unpack
from osgeo import gdal
from osgeo import osr
 
# inputfile
BinaryFile = "/path/to/raw6hr0006100.bin" # '3B42_daily.2009.05.31.7.bin'
 
# open binary file
f = open(BinaryFile, "rb")
 
# set file dimensions
xs = 1440
ys = 400
 
# set number of bytes in file
NumbytesFile = xs * ys
 
# number of columns in row
NumElementxRecord = -xs
 
# create empty array to put data in
myarr = []
 
# loop trough the binary file row by row
for PositionByte in range(NumbytesFile,0, NumElementxRecord):
 
        Record = ''
 
        # the dataset starts at 0 degrees, use 720 to convert to -180 degrees
        for c in range (PositionByte-720, PositionByte, 1):
                f.seek(c * 4)
                DataElement = unpack('>f', f.read(4))
                Record = Record  + str("%.2f" % DataElement + ' ')
 
        # 0 - 180 degrees
        for c in range (PositionByte-1440 , PositionByte-720, 1):
                f.seek(c * 4)
                DataElement = unpack('>f', f.read(4))
                Record = Record  + str("%.2f" % DataElement + ' ')
 
        # add data to array
        myarr.append(Record[:-1].split(" "))
 
# close binary file
f.close()
 
# Array to numpy float
myarr = np.array(myarr).astype('float')
 
# set values < 0 to nodata
myarr[myarr < 0] = -9999
 
# mirror array
myarr = myarr[::-1]
 
# define output name
outname = "/path/to/persian.tif"
 
# set coordinates
originy = 50
originx  = -180
pixelsize = 0.25
transform= (originx, pixelsize, 0.0, originy, 0.0, -pixelsize)
driver = gdal.GetDriverByName( 'GTiff' )
 
# set projection
target = osr.SpatialReference()
target.ImportFromEPSG(4326)
 
## write dataset to disk
outputDataset = driver.Create(outname, xs,ys, 1,gdal.GDT_Float32)
outputDataset.SetGeoTransform(transform)
outputDataset.SetProjection(target.ExportToWkt())
outputDataset.GetRasterBand(1).WriteArray(myarr)
outputDataset.GetRasterBand(1).SetNoDataValue(-9999)
outputDataset = None