一、前言
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
