用MATLAB/Python/hdp 打开TRMM HDF格式文件

详情在README.TRMM_V7.pdf的4.1章节中。

MATLAB

需要安装额外的Mapping Toolbox
MATLAB can be used to load, manipulate, and view TRMM precipitation data. To load the
precipitation variable from the aforementioned TRMM file into MATLAB type:

precip = permute(hdfread('3B43.19990101.7.HDF','precipitation'),[2 1]);

This will load the data into a matrix called precip. Missing data are represented by -9999.9, but MATLAB doesn’t know that this value refers to missing data. The simplest way to replace the missing numeric values with MATLAB’s not-a-number (NaN) values, is to type:

precip(precip < 0) = NaN;

figure;

axesm('MapProjection','eqdcylin','maplatlimit',[-50 50],'maplonlimit',[-180 180],'ParallelLabel','on','PlabelMeridian','west','MeridianLabel','on','MLabelParallel','south','FontSize',6,'FontWeight','bold','PLineLocation',20,'MLineLocation',20);

latitudes = -49.875:0.25:49.875;
longitudes = -179.875:0.25:179.875;

[latGrid, lonGrid] = meshgrat(latitudes,longitudes);
geoshow(latGrid, lonGrid, double(precip), 'DisplayType','texturemap');
caxis([0 5]);
colormap(flipud(hot(21)));

% This line places the colorbar
chandle = colorbar('Location','EastOutside','FontSize',6,'FontWeight','bold'); 
 % Set the colorbar's label
set(get(chandle,'ylabel'),'String','Average Rain Rate (mm/h)','FontSize',10,'FontWeight','Bold');
set(chandle,'YTick',0:5);

states = geoshape(shaperead('landareas', 'UseGeoCoords', true)); 
geoshow(states,'DefaultFaceColor','none','DefaultEdgeColor','k'); 
tightmap
title('July 1999 Rain Rate (mm/h)','FontSize',8,'FontWeight','bold'); 

Python

python 需要 numpy, matplotlib, basemap, 和 pyhdf, 其中可以参考basemap的库安装

# This is a test script that reads and plots the TRMM 3B42 3-hourly data. from mpl_toolkits.basemap import Basemap, cm
import matplotlib.pyplot as plt
import numpy as np
from pyhdf.SD import SD, SDC
dataset = SD('/path/to/3B42.20120824.12.7.HDF', SDC.READ)
precip = dataset.select('precipitation') precip = precip[:]
precip = np.transpose(precip)
theLats = np.arange(-49.875,50,0.25) theLons = np.arange(-179.875,180,0.25)
# Set all the missing values less than 0 to NaNs np.putmask(precip,precip<0,np.nan)
# Plot the figure, define the geographic bounds fig = plt.figure(dpi=300)
latcorners = ([-50,50])
loncorners = ([-180,180])
m = Basemap(projection='cyl',\ llcrnrlat=latcorners[0],urcrnrlat=latcorners[1],llcrnrlon=loncorners[0],urcrnrlon=loncorners[1])
# Draw coastlines, state and country boundaries, edge of map. m.drawcoastlines()
m.drawstates()
m.drawcountries()
# Draw filled contours.
clevs = np.arange(0,5.01,0.5)
# Define the latitude and longitude data
x, y = np.float32(np.meshgrid(theLons, theLats))
cs = m.contourf(x,y,precip,clevs,cmap=cm.GMT_drywet,latlon=True)
parallels = np.arange(-50.,51,25.) m.drawparallels(parallels,labels=[True,False,True,False]) meridians = np.arange(-180.,180.,60.) m.drawmeridians(meridians,labels=[False,False,False,True])
     64
 # Set the title and fonts
plt.title('24 August 2012 1200 UTC Rain Rate')
font = {'family' : 'normal', 'weight' : 'bold', 'size' : 6} plt.rc('font', **font)
# Add colorbar
cbar = m.colorbar(cs,location='right',pad="5%") cbar.set_label('mm/h')
plt.savefig('testTRMMmap.png',dpi=300)

hdp and ncdump

The HDF Toolkit ships with two binary executables, hdp and ncdump, that can be used to extract values from any HDF file. These are also available as standalone executable from the utilities folders found within each operating system at:
ftp://ftp.hdfgroup.org/HDF/HDF_Current/bin.

ncdump can only read HDF files if your local copy of netCDF was originally compiled with HDF support.
To dump the entire file: hdp <filename> or ncdump <filename>
To get just the header information: hdp dumpsds -h <filename> or ncdump -h <filename>
A partial example of output from hdp dumpsds -h 3B42.20120824.12.7.HDF is given below. (The ncdump -h output is similar.)

File attributes:
Attr0: Name = FileHeader
Type = 8-bit signed char
Count= 357
Value = AlgorithmID=3B42;\012AlgorithmVersion=3B4
2_7.0;\012FileName=3B42.20120824.12.7.HDF ;\012GenerationDateTime=2012-10-26T14:07: 33.000Z;\012StartGranuleDateTime=2012-08- 24T10:30:00.000Z;\012StopGranuleDateTime= 2012-08-24T13:29:59.999Z;\012GranuleNumbe r=;\012NumberOfSwaths=0;\012NumberOfGrids =1;\012GranuleStart=;\012TimeInterval=3_H OUR;\012ProcessingSystem=PPS;\012ProductV ersion=7;\012MissingData=;\012
Attr1: Name = FileInfo
Type = 8-bit signed char
Count= 253
Value = DataFormatVersion=m;\012TKCodeBuildVersio
n=1;\012MetadataVersion=m;\012FormatPacka ge=HDF Version 4.2 Release 4, January 25,
2009;\012BlueprintFilename=TRMM.V7.3B42. blueprint.xml;\012BlueprintVersion=BV_13; \012TKIOVersion=1.6;\012MetadataStyle=PVL ;\012EndianType=LITTLE_ENDIAN;\012
Attr2: Name = GridHeader Type = 8-bit signed char
Count= 231
Value = BinMethod=ARITHMETIC_MEAN;\012Registratio
n=CENTER;\012LatitudeResolution=0.25;\012 LongitudeResolution=0.25;\012NorthBoundin gCoordinate=50;\012SouthBoundingCoordinat e=-50;\012EastBoundingCoordinate=180;\012 WestBoundingCoordinate=-180;\012Origin=SO UTHWEST;\012
Variable Name = precipitation Index = 0
Type= 32-bit floating point Ref. = 2
Compression method = NONE Rank = 2
Number of attributes = 1 Dim0: Name=nlon
Size = 1440
Scale Type = number-type not set Number of attributes = 0
Dim1: Name=nlat Size = 400
Scale Type = number-type not set
Number of attributes = 0 Attr0: Name = units
Type = 8-bit signed char
 66
 Count= 5 Value = mm/hr
... and so on ... This will list all of the variables in the same manner.