Hi Everyone,
I've written a Python script to read a binary TRMM data file and write
the data, with coordinates, to a CSV file so that I can import the
data into a GIS.
However, I'd like some doublde check on whether the script assigns the
correct precipitation value to each coordinate pair.
The NASA documentation states:
"It is size 1440x480, with X (longitude) incrementing most rapidly
West to East from the Prime
Meridian, and then Y (latitude) incrementing North to South from the
northern edge."
The Python code is (thanks to Scott Sinclair for his help with this):
[code]
# -*- coding: utf-8 -*-
'''
Created on Tue Apr 20 15:39:57 2010
Read TRMM data from a binary file that covers the earth from 60°N
to 60°S, construct the lat-long coordinates for each cell and write
the
coordinates and precipitation for each cell to a CSV file.
TRMM data are in 0.25 degree squares, resulting in 480 rows and
1440 columns of data.
'''
import os
import numpy as np
precip_scale_factor = 100.0
rows = 480
cols = 1440
def read_trmm_data(fname):
''' read the data from the file 'fname' to a 480x1440 array and
return the array
'''
fp = open(fname, 'rb')
data_string = fp.read()
fp.close()
precip = np.fromstring(data_string[2880:1385280], np.int16)
precip = precip.byteswap()
precip = np.asarray(precip, np.float32)
precip /= precip_scale_factor
precip = precip.reshape(rows, cols)
def write_trmm_data():
''' build 2 480x1440 arrays for the cell coordinates, one
for lattitude and one for longitude and combine the elements
from these two arrays with the elements in the
trmm_precipitation
array to write a CSV file with the elements in
the order lat, long, precipitation
'''
# create one rwo of N-S coordinates
# coordinate locates the centre of the cell
lat=np.arange(59.875,-60.125,-0.25, dtype=float)
z=np.arange(59.875,-60.125,-0.25, dtype=float)
for i in range(1,1440):
lat=np.vstack((lat,z))
lat=lat.transpose()
# create one row of E-W coordinates
# coordinate locates the centre of the cell
lng=np.arange(0.125,360.125,+0.25, dtype=float)
z=np.arange(0.125,360.125,+0.25, dtype=float)
for i in range(1,480):
lng=np.vstack((lng,z))
#open file and write the data
out_text = file(fname+".txt","w")
headings='lat,long,precipitation\n'
out_text.write(headings)
for i in np.ndindex(precip.shape):
filerow='%g,%g,%g\n' % (lat[i],lng[i],precip[i])
out_text.write(filerow)
out_text.close()
fname="3B42RT.2010032409.6.bin"
read_trmm_data(fname)
write_trmm_data()
[/code]
This code gives an output file with the entries (the -319.99 indicates
no data):
[output]
lat,long,precipitation
59.875,0.125,-319.99
59.875,0.375,-319.99
59.875,0.625,-319.99
59.875,0.875,-319.99
59.875,1.125,-319.99
.
.
.
-59.875,358.875,-319.99
-59.875,359.125,-319.99
-59.875,359.375,-319.99
-59.875,359.625,-319.99
-59.875,359.875,-319.99
[/output]
I can't think of a way to test this output. Can someone perhaps help
me?
Thanks
Hanlie
--
Subscription settings:
http://groups.google.com/group/python-gis-sig/subscribe?hl=en
