Matlab读写HDF-EOS文件

2010年11月18日星期四 0 评论
最近才发现Matlab集成了读写HDF文件(包括HDF-EOS文件)的功能。

HDF-EOS的三种数据类型及其对应的Matlab函数如下:
grid(GD) - hdfgd
point(PT) - hdfpt
swath(SW) - hdfsw

此外,Matlab还提供了GUI的HDF文件导入工具,可以通过 “hdftool”命令启动。

hdfgd的帮助需要用 “help hdfgd" 查看。
hdfgd可以认为是HDF-EOS库中所有Grid函数的一个共用接口,即通过hdfgd来调用不同的Grid函数。
以下是读取一个HDF-EOS文件的代码示例,可以看出完成相同的任务,Matlab比C/C++所需代码量少得多,当然执行效率会不如C/C++。

hdf = 'MOD15A2.A2006361.h10v04.005.2008134052818.hdf';
fileid = hdfgd('open', hdf, 'DFACC_READ');                 %对应GDopen
gridid = hdfgd('attach', fileid, 'MOD_Grid_MOD15A2');      %对应GDattach
[data, status] = hdfgd('readfield', gridid, 'Lai_1km', [], [], []);  %对应GDreadfield
hdfgd('detach', gridid);  %对应GDdetach
hdfgd('close', fileid);   %对应GDclose

全球LAI产品(Global LAI Products)

2010年11月14日星期日 0 评论
文献[1]对4种覆盖全球的LAI产品进行了综合评价MODISCYCLOPESECOCLIMAPGLOBCARBON,如下表所示。其它还有一些全球LAI产品,如POLDERLand-SAF。还有一些区域性的LAI产品,如加拿大范围的LAI产品。这些产品中,既有基于物理模型,实现时采用查找表(LUTLook-Up Table)或神经网络;也有基于LAI与植被指数(VIVegetation Index)统计关系的。

LAI产品名
时间分辨率
空间分辨率
年份
传感器
产品生成算法
MODIS
8
1km
2000-
MODIS
主算法:3维辐射传输模型(LUT实现)
备用算法:基于NDVI的经验关系
CYCLOPES
10
1km
1999-2007
VEGETATION
1维辐射传输模型(神经网络实现)
ECOCLIMAP
1
1km
climatology
AVHRR
基于NDVI的经验关系
GLOBCARBON
1
0.25°, 0.5°, 10km
1998-2007
AATSR, ATSR, MERIS, VEGETATION
基于VI的统计关系

(climatology的意思是:the mean state based on a historical set of observations)

以上四种产品的相互比较及与地面实测数据的比较[1]表明,MODIS和CYCLOPES对真实地表及其空间差异的表达,LAI时间序列特征及年份间的差异都优于ECOCLIMAP和GLOBCARBON。但另一方面,MODIS和CYCLOPES的反演成功率却低于ECOCLIMAP和GLOBCARBON,因此MODIS和CYCLOPES存在较多缺失值,时间序列不够连续,影响了产品的应用。

事实上,从上表也可以看出,目前只有MODIS的LAI产品在持续发布,其它产品处于产品评价和算法改进过程中。

使用HDF库读写HDF4文件(C/C++)

2010年8月27日星期五 0 评论
一、引言

Hierarchical Data Format (HDF)是由HDF Group开发维护的一种开放文件格式,目前HDF有HDF4和HDF5两个版本,它们相互并不兼容。相应的HDF-EOS也有两个版本,基于HDF4的叫做HDF-EOS2,基于HDF5的叫做HDF-EOS5,MODIS数据产品使用的是HDF-EOS2文件格式。HDF Group提供了C/C++,Fortran,Java库以对HDF文件进行读写操作HDF4最新版本号是4.2.5,并针对多种平台进行了预先编译。本文讨论如何用官方编译好的Windows平台的C/C++库中的SD API对HDF4文件进行读写操作。

本文参考以下三个文档,均可从HDF Group网站获得
HDF4 User’s Guide
HDF4 Reference Manual
sd_rd.c

二、HDF文件格式

图片引自《HDF4 User’s Guide》
如上图所示,一个HDF文件中可以有6种数据结构:Raster Image,Palette,Scientific Data Set,Annotation,Vdata和Vgroup。与这六种数据结构相应的有5类API:

SD API - 对Scientific Data Set进行读写操作
VS API - 对Vdatas进行读写操作
V API - 对Vgroups进行读写操作
GR API - 对Raster Images和Palettes进行读写操作
AN API Stores - 对Annotations进行读写操作

说到这里则有必要提及两个现成的HDF查看软件:HDF Explorer和HDFView,HDF Explorer是私有软件,而HDFView是由HDF Group开发维护的。HDF Explorer默认支持HDF-EOS,HDFView则需要安装HDF-EOS_Plugin才能获得较好的HDF-EOS支持。下面是HDFView截图。


在HDFView中查看HDF文件中DataFields的Metadata
红框标示了数据集的类型和数据的类型

三、例子
本例子读入一个MOD15A2 HDF文件的LAI数据(uint8),将其除以10(float)再输出为新的HDF文件。在程序中使用数组来存储原始LAI数据和除以10的LAI数据,数组类型需要与相应HDF文件中指定的数据类型严格对应,否则会发生数据读写错误。

完整的程序(ReadWrite_SDS.c),makefile,数据文件MOD15A2.A2006361.h10v04.005.2008134052818.hdf。需要链接的库文件可以查阅下载下来的HDF库中的INSTALL_WINDOWS.txt。在本例子程序中,只链接了hd425m.lib和hm425m.lib,这可以从makefile中看出来。

读入HDF文件的代码如下:
/* Read data to LAI_MODIS array */
sd_id = SDstart("MOD15A2.A2006361.h10v04.005.2008134052818.hdf", DFACC_READ);
sds_index = SDnametoindex(sd_id, "Lai_1km");
sds_id = SDselect (sd_id, sds_index);

输出HDF文件代码如下:

/* Create and open the file and initiate the SD interface. */
sd_id = SDstart("MyMOD15A2.hdf", DFACC_CREATE);
sds_id = SDcreate(sd_id, "LAI_divided", DFNT_FLOAT32, rank, dims);
istat = SDwritedata(sds_id, start, NULL, edges, (VOIDP)LAI);

MODIS产品笔记整理

2010年8月2日星期一 0 评论
之前将笔记零碎的记录在Google Docs里,现在把它整理一下。

一、关于MODIS命名

我们知道,Terra和Aqua卫星都搭载有MODIS传感器,在MODIS产品名和文件名中则反映了数据来源。

MOD开头表示是用Terra卫星数据生成的产品,
MYD开头表示是用Aqua卫星数据生成的产品,
MCD开头则表示是结合了两颗卫星数据合成的产品,字母C表示Combined。

完整的MODIS陆面产品列表参考以下网址:

二、关于数据格式

一般使用的MODIS数据产品有两种:

.hdf文件:其文件扩展名为.hdf,但严格的说,它的文件格式是HDF-EOS。HDF-EOS是HDF的扩展,包含了geo-location,投影信息等。以.hdf文件提供MODIS的数据产品,每个文件就是经过Sinusoidal Tiling System分片后的一个tile。一个tile覆盖1200km*1200km。.hdf文件的陆面数据产品可以从https://lpdaac.usgs.gov获得。

.txt文件:MODIS数据产品的子集,覆盖Selected Sites(一般是通量站点)7km*7km的大小。.hdf文件是根据产品的发布时间,如MOD09A1每个tile每8天提供一个HDF文件,而.txt则是一个站点一个产品一个文件,因此.txt文件可以更方便的得到时间序列。.txt文件的陆面数据产品可以从http://daac.ornl.gov获得。

三、 关于QC或是QA

MODIS每种数据产品都提供了质量控制(QC: Quality Control)或质量保证(QA: Quality Assurance)字段。QC信息是像元级的,标识了数据质量的高低,用户根据实际需求,使用QC信息筛选出质量足够好的数据。下面分别对MOD15A2和MOD09A1产品的QC字段进行说明。

• MOD15A2产品有两个QC字段:FparLai_QC和FparExtra_QC,FparLai_QC(8位)的描述如下:

Bit No.
Parameter Name
Bit Comb.
FparLai_QC
0MODLAND_QC bits0Good quality (main algorithm with or without saturation)
1Other Quality (back-up algorithm or fill values)
1Sensor0Terra
1Aqua
2DeadDetector0Detectors apparently fine for up to 50% of channels 1,2
1Dead detectors caused >50% adjacent detector retrieval
3–4CloudState (inherited from Aggregate_QC bits {0,1} cloud state)000 Significant clouds NOT present (clear)
011 Significant clouds WERE present
102 Mixed cloud present on pixel
113 Cloud state not defined, assumed clear
5–7SCF_QC (five level confidence score)0000, Main (RT) method used, best result possible (no saturation)
0011, Main (RT) method used with saturation. Good,very usable
0102, Main (RT) method failed due to bad geometry, empirical algorithm used
0113, Main (RT) method failed due to problems other than geometry, empirical algorithm used
1004, Pixel not produced at all, value coudn't be retrieved (possible reasons: bad L1B data, unusable MODAGAGG data)
(此表引用自:https://lpdaac.usgs.gov/lpdaac/products/modis_products_table/leaf_area_index_fraction_of_photosynthetically_active_radiation/8_day_l4_global_1km/mod15a2

MOD15 ATBD中指出:The key indicator of retrieval quality of the LAI/FPAR product is SCF_QC bitfield.

MOD09A1产品也包含了两个QC字段:500m Reflectance Band Quality和500m State Flags,500m Reflectance Band Quality(32位)的描述如下:


Bit No.
Parameter Name
Bit Comb.
Sur_refl_qc_500m
31adjacency correction performed1yes
0no
30atmospheric correction performed1yes
0no
26–29band 7 data quality four bit range0000highest quality
1000dead detector; data interpolated in L1B
1001solar zenith >= 86 degrees
1010solar zenith >= 85 and <>
1011missing input
1100internal constant used in place of climatological data for at least one atmospheric constant
1101correction out of bounds pixel constrained to extreme allowable value
1110L1B data faulty
1111not processed due to deep ocean or clouds
22–25band 6 data quality four bit rangeSAME AS BAND ABOVE
18–21band 5 data quality four bit rangeSAME AS BAND ABOVE
14–17band 4 data quality four bit rangeSAME AS BAND ABOVE
10–13band 3 data quality four bit rangeSAME AS BAND ABOVE
6–9band 2 data quality four bit rangeSAME AS BAND ABOVE
2–5band 1 data quality four bit rangeSAME AS BAND ABOVE
0–1MODLAND QA bits00corrected product produced at ideal quality all bands
01corrected product produced at less than ideal quality some or all bands
10corrected product not produced due to cloud effects all bands
11corrected product not produced due to other reasons some or all bands may be fill value [Note that a value of (11) overrides a value of (01)].
(此表引用自:https://lpdaac.usgs.gov/lpdaac/products/modis_products_table/surface_reflectance/8_day_l3_global_500m/mod09a1

在我的工作中,只使用了MODLAND QA bits来进行质量控制。

这里需要特别指出的是:在QC字段中第0位是低位,即最右边那一位。

四、 地表覆盖类型(Land Cover Type)数据产品MCD12Q1

MCD12Q1是结合了Terra和Aqua两颗卫星的数据生成的地表覆盖类型数据产品,时间分辨率1年,空间分辨率500m。
MCD12Q1中包含5种分类体系(classification schemes)的分类结果,5种分类体系如下:

  • Land Cover Type 1: IGBP global vegetation classification scheme
  • Land Cover Type 2: University of Maryland (UMD) scheme
  • Land Cover Type 3: MODIS-derived LAI/fPAR scheme
  • Land Cover Type 4: MODIS-derived Net Primary Production (NPP) scheme
  • Land Cover Type 5: Plant Functional Type (PFT) scheme


ClassIGBP (Type 1)UMD (Type 2)LAI/fPAR (Type 3)NPP (Type 4)
0WaterWaterWaterWater
1Evergreen Needleleaf forestEvergreen Needleleaf forestGrasses/Cereal cropsEvergreen Needleleaf vegetation
2Evergreen Broadleaf forestEvergreen Broadleaf forestShrubsEvergreen Broadleaf vegetation
3Deciduous Needleleaf forestDeciduous Needleleaf forestBroadleaf cropsDeciduous Needleleaf vegetation
4Deciduous Broadleaf forestDeciduous Broadleaf forestSavannaDeciduous Broadleaf vegetation
5Mixed forestMixed forestEvergreen Broadleaf forestAnnual Broadleaf vegetation
6Closed shrublandsClosed shrublandsDeciduous Broadleaf forestAnnual grass vegetation
7Open shrublandsOpen shrublandsEvergreen Needleleaf forestNon-vegetated land
8Woody savannasWoody savannasDeciduous Needleleaf forestUrban
9SavannasSavannasNon-vegetated
10GrasslandsGrasslandsUrban
11Permanent wetlands
12CroplandsCroplands
13Urban and built-upUrban and built-up
14Cropland/Natural vegetation mosaic
15Snow and ice
16Barren or sparsely vegetatedBarren or sparsely vegetated
254UnclassifiedUnclassifiedUnclassifiedUnclassified
255Fill ValueFill ValueFill ValueFill Value
(此表引用自:https://lpdaac.usgs.gov/lpdaac/products/modis_products_table/land_cover/yearly_l3_global_500_m/mcd12q1