Working with TRMM rainfall data

In the next exercise, you are going to download and extract TRMM data.

This data can be obtained in HDF and NETCDF format. We will download data in the NETCDF format.

In this assignment you will create some python scripts to open the files and extract the data. For this training, you will just execute the scripts and save the data. It is not necessary yet to understand all regarding the scripts yet.

step 1: Go to this website and download the monthly TRMM data in NETCDF format.

Question 4.1: What is the name and temporal resolution of the product?

step 2: Save the files in d:\training1\trmm\

step 3: Open the Osgeo4W shell.

step 4: Go to the directory:

d:
cd training1\trmm\

in the next step you are going to work with the Geospatial Data Abstraction Library.

Geospatial Data Abstraction Library (GDAL/OGR) is a cross platform C++ translator library for raster and vector geospatial data formats that is released under an X/MIT style Open Source license by the Open Source Geospatial Foundation. As a library, it presents a single abstract data model to the calling application for all supported formats. It also comes with a variety of useful commandline utilities for data translation and processing.

GDAL supports over 50 raster formats, and OGR over 20 vector formats.

It provides the primary data access engine for many applications including MapServer, GRASS, QGIS, and OpenEV. It is also utilized by packages such as OSSIM, Cadcorp SIS, FME, Google Earth, VTP, Thuban, ILWIS, MapGuide and ArcGIS. GDAL/OGR is the most widely used geospatial data access library [ref].

step 5: Get the information of the netcdf file by typing

 gdalinfo 3B43.20000101.7A.nc 

gdalinfo returns you the information on datalayer. But GDAL also includes numerous other functions.

Question 4.2: How many data layers does the file have?

Now you are going to use Geany. Geany is a lightweight GUI text editor, including basic Integrated Development Environment (IDE) features. It is designed to have short load times.

step 6: Open Geany and open a new file.

The easiest way to handle netCDF data is to open the data in script such as Python and then extract the relevant information. In the next steps you are going to create such a script.

step 7. Save the file as script4a.py in d:\training1\scripts\

step 8: Add the following text to the script:

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

# location to the inputfile
fn = "/path/to/vietnamtraining/trmm/3B43.20000101.7A.nc"

# open the file
dat = gdal.Open(fn, gdal.GA_ReadOnly)

# Get the Precipitation dataset
dataset = dat.GetSubDatasets()[0]
print dataset

The first section opens all libraries. a library is a collection of precompiled routines that a program can use.

Next the location of the file is specified.

step 9: Adjust the name of the file to the right input location.

Then the file is opened using the GDAL libary.

In the last step the name of the dataset is printed.

Step 10. Go to the Osgeo4W shell and type

 cd d:\training1\scripts\

Step 11: run the script by typing

 python script4a.py 

question 4.3: What is the output? What does pcp stand for? what is the data format?

The variable precip contains the two dimensional matrix of the precipitation data. With print precip.GetMetadata() you can print the metadata.

step 12: create a new script4b.py, paste the code below in the script and execute the code.

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

# location to the inputfile
fn = "/path/to/3B43.20000101.7A.nc"

# open the file
dat = gdal.Open(fn, gdal.GA_ReadOnly)

# Get the Precipitation dataset
dataset = dat.GetSubDatasets()[0]

# read the precipitation dataset
precip = gdal.Open(dataset[0],gdal.GA_ReadOnly)

# print metadata
print precip.GetMetadata()

question 4.4: What is metadata?

question 4.5. What is the unit of the data?

question 4.6: What is the nodata (fill) value?

step 13: Create a new script called script4c.py and display the data by executing the following script:

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

# location to the inputfile
fn = "/path/to/3B43.20000101.7A.nc"

# open the file
dat = gdal.Open(fn, gdal.GA_ReadOnly)

# Get the Precipitation dataset
dataset = dat.GetSubDatasets()[0]

# read the precipitation dataset
precip = gdal.Open(dataset[0],gdal.GA_ReadOnly)

# get the data of the precipitation dataset
preciparr = precip.ReadAsArray()

# show the data
plt.imshow(preciparr)
plt.colorbar()
plt.show()

Here it is important to note that the data is stored as a matrix (or array) in the variable preciparr.

A one dimenstions array can looks like this: myArray=[1,2,3,4,5,6]

A two dimensional array can looks like this: myArray=[[1,2],
[3,4],
[5,6]]

As such, data can be displayed in and image e.g.: [[12,35,11,12],
[32,55,11,59],
[22,25,11,12],
[24,15,91,72]]

Question 4.7: Save the image and add it to the report.

The script below writes the precipitation data to a Geotiff file on your disk.

step 14: Create a new script4d.py and execute the text below. Adjust the output path.

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

# location to the inputfile
fn = "/path/to/3B43.20000101.7A.nc"

# open the file
dat = gdal.Open(fn, gdal.GA_ReadOnly)

# Get the Precipitation dataset
dataset = dat.GetSubDatasets()[0]

# read the precipitation dataset
precip = gdal.Open(dataset[0],gdal.GA_ReadOnly)

# get the data of the precipitation dataset
preciparr = precip.ReadAsArray()

# get geotransform
transform = precip.GetGeoTransform()

# set geotif driver
driver = gdal.GetDriverByName( 'GTiff' )

# get x,y dimensions of the map
xs = precip.RasterXSize
ys = precip.RasterYSize

# set output name
outname = "/path/to/january.tif"

# 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(preciparr)
outputDataset.GetRasterBand(1).SetNoDataValue(-9999)
outputDataset = None

step 15: Open the Geotiff file in Qgis.

Question 4.7: What is the range of values?

Question 4.8: What should be done to convert the unit to mm/month?

step 16: Use the script to extract all files for the year 2000.

step 17:. Use the Vietnam shapefile to create a monthly precipitation map for Vietnam for each month.

Question 4.9: make a screenshot of the result (December) and include it in the report.

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

Month	Min	Mean	Max
Jan
Feb
..		.

Question 4.10: include the table in the report.

Question 4.11: Compare the averaged yearly precipitation rates of CHIRPS with TRMM, what is the difference?