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?
What should be done to convert the TRMM 3B43 mm/hr unit to mm/month?
LikeLike
By summing it for a month
LikeLike
Could you provide a worked example for converting the rows and columns of original TRMM 3b43 data?Thank you very much.
LikeLike