Creating datasets with monthly rainfall quantities

From 6-hourly to monthly rainfall datasets using python

In this exercise you worked with 6-hourly datasets. However, this resolution is not necessarily required for many applications. In this exercise, you will combine all 6-hourly maps into one monthly map.

We can do this by
1. Listing all files in a folder.
2. Create a map with only zeros.
3. Opening the files one by one.
4. Adding the opened maps to the map with zeros.

This can be done using a loop.

step 1. Create script6a.py with the following code:

# import libraries
import os

# input directory
directory = "/path/to/persiann/"

# list files in directory
filesindir = os.listdir(directory)

# print files in folder
for files in filesindir:
	print files

Question 6.1. Describe in your own words what happens in this script.

step 2: We only want the binary files and not the compressed files. Create script6b.py and use the code below:

# import libraries
import os

# input directory
directory = "/path/to/persiann/"

# list files in directory
filesindir = os.listdir(directory)

# print files in folder
for files in filesindir:
	# try except statement to catch errors
	try:
		# split the filename based on the .
		splitfilename = files.split(".")
		# if the length of the name is longer than 3 [name, tar, gz] then pass else print
		if len(splitfilename) == 2:
			print files
	except:
		pass

question 6.2. Compare the output with the output of the previous script: what has changed?

step 3. Now we will combine the script above with the script of the previous exercise. The new script (script6c.py) will combine all monthly maps into a monthly product. Create the new script and run the code below.

# import libraries
import matplotlib.pyplot as plt
import numpy as np
from struct import unpack
from osgeo import gdal
from osgeo import osr
import os

def ReadBinaryData(fname):

	# open binary file
	f = open(directory + fname, "rb")

	# 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()
	myarr = np.array(myarr).astype('float')

	return myarr

# set file dimensions
xs = 1440
ys = 400

# set number of bytes in file
NumbytesFile = xs * ys

# number of columns in row
NumElementxRecord = -xs

# create array to put data in
monthlydata = np.zeros(shape=(ys,xs))

# input directory
directory = "/path/to/persiann/"

# list files in directory
filesindir = os.listdir(directory)

# print files in folder
for files in filesindir:
	# try except statement to catch errors
	try:
		# split the filename based on the .
		splitfilename = files.split(".")
		# if the length of the name is longer than 3 [name, tar, gz] then pass else print
		if len(splitfilename) == 2 and splitfilename[1] == "bin":
			print files
                        data = ReadBinaryData(files)
                        data[data < 0] = 0
			monthlydata += data
	except:
		pass

# set values < 0 to nodata
monthlydata[monthlydata < 0] = -9999

# mirror array
monthlydata = monthlydata[::-1]

# define output name
outname = "/path/to/March2000.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(monthlydata)
outputDataset.GetRasterBand(1).SetNoDataValue(-9999)
outputDataset = None

Question 6.3: Explain again in your own words what the script does.

step 4: Clip Vietnam from the newly created map.

<!–
step 5: Do the same for April – Dec 2000.

step 6: right click the maps, go to properties and metadata. Use the statistics data to fill out the following table for each month: