Working with PERSIANN rainfall data

Processing near real-time rainfall information on precipitation

PERSIANN, “Precipitation Estimation from Remotely Sensed Information using Artificial Neural Networks”,is a satellite-based precipitation retrieval algorithm that provides near real-time rainfall information.

The algorithm uses infrared (IR) satellite data from global geosynchronous satellites as the primary source of precipitation information. Precipitation from IR images is based on statistical relationship between cloud top temperature and precipitation rates. The IR-based precipitation estimates are then calibrated using satellite microwave data available from low Earth orbit satellites (e.g., Tropical Rainfall Measuring Mission Microwave Imager, Special Sensor Microwave Imager, Advanced Microwave Scanning Radiometer‐Earth observing system). The calibration technique relies on an adaptive training algorithm that updates the retrieval parameters when microwave observations become available (approximately at 3 hours intervals). [1]

The PERSIANN data can be acquired in Binary format. This format is different from the geotif, netCDF and HDF we have seen before. As such, a different approach is need to extract the data.

step 1: Go to the ftp site and download the data of March 2000.

step 2: Extract and decompress the files in the folder d:/persiann/

When extracted using tar and then decompressed , the 6-hrly files will have the following names: raw6hrYYDOYHS.bin

Where:
YY : 2 digit year starting from 00 for 2000
DOY: Day of year
HS: Starting hour of the 6 hours accumulation period, (0,6,12,18)
for example, the file raw6hr0512518.bin represents the precipitation accumulated between 18 and 24, on May 5th, 2005

The data is stored in C style row centric format. with the first value being centered at 49.875, 0.125, the second value at 49.875, 0.375, and the last 2 values are centered at: -49.875, 359.625 & -49.875, 359.875

Data are in 4-byte binary float from a SUN system (big-endian).

The following step is required to extract the binary data from the binary file. It basically scans the data row by row and puts the data in a matrix.

step 3:. Create the script script5a.py in d:\training1\scripts\ and run the code below.

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

If the image with rainfall appears you have successfully imported the binary data into a numpy array.

Question 5.1: Save the image to a file and include it in the report.

Question 5.2: Explain in 5 simple lines what the script does.

4. Create the script script5b.py in d:\training1\scripts\ and run the code below. Make sure you correct the input and output path.

# 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

Question 5.3: Review the script and explain in your own words how the script works.

Question 5.4. Open the file you just created and answer the following questions:
– cell size:
– extent:
– maximum value:
– mean value:
– minimum value:
– Dimensions

Question 5.5: Now you have extracted one binary file into a GeoTIFF. What should be done to create a monthly data product?