MATLAB处理GPM(IMERG/GSMaP)卫星降水数据
2019-04-29 本文已影响0人
zengsk
概述
GPM卫星群(来源于NASA官网)
GPM 介绍
2014年2月28日发射的全球降水观测计划（GPM）标志着卫星降水正式由TRMM时代跨入了GPM时代。GPM主卫星携带了最新的Ku/Ka多普勒双频降雨雷达（DPR）和微波成像仪（GMI），将为遥感水文科学社区提供更高时空分辨率（1-hour, 10-km）和更高精度的全球卫星降水观测。GPM计划主要由美国宇航局（NASA）和日本宇宙航天机构（JAXA）联合实施，因此GPM时代主推的两套新一代的卫星降水便是美国宇航局负责研发的IMERG和日本宇宙航天机构负责的GSMaP。
GPM-IMERG
- 数据介绍
- IMERG是专为全球降水计划GPM而生的最新一代多卫星融合反演降水数据，是GPM的3级产品。它充分利用GPM平台上所有的卫星传感器提供的数据（包括主被动微波传感器和各类红外数据传感器等等），也充分借鉴之前TRMM时代基本成熟的各类卫星降水反演算法进行有机融合。
-
IMERG目前提供三套类型的卫星降水数据，分别是Early，Late，Final三个版本。IMERG生成系统在实时阶段运行一次后得到Early产品，再运行一次后得到Late产品，在这过程中，最主要的不同在于Early中只采用了云移动矢量传播算法中的前向传播算法，而Late在此基础上还加用了后向传播算法。在滞时产品Final的处理中，在Late的基础上引进了更多的传感器数据源，包括引入了全球雨量站点进行校正。
GPM IMERG降水数据介绍
-
数据下载方法(ftp下载)
作者常用的是FlashFXP、cuteFTP、ftpdownload软件进行批量下载！！
ftp下载数据
-
GPM数据的组织结构及其变量
可以用matlab的h5disp(filename) 或者 h5info("dataset")进行查看
% HDF5 3B-HHR.MS.MRG.3IMERG.20170601-S003000-E005959.0030.V05B.HDF5
% '
% 'FileInfo': 'DataFormatVersion=cn;
% TKCodeBuildVersion=3;
% MetadataVersion=cv;
% FormatPackage=HDF5-1.8.9;
% BlueprintFilename=GPM.V1.3IMERGHH.blueprint.xml;
% BlueprintVersion=BV_58;
% TKIOVersion=3.80.33;
% MetadataStyle=PVL;
% EndianType=LITTLE_ENDIAN;
% Group '/Grid'
% Attributes:
% 'GridHeader': 'BinMethod=ARITHMETIC_MEAN;
% Registration=CENTER;
% LatitudeResolution=0.1;
% LongitudeResolution=0.1;
% NorthBoundingCoordinate=90;
% SouthBoundingCoordinate=-90;
% EastBoundingCoordinate=180;
% WestBoundingCoordinate=-180;
% Origin=SOUTHWEST;
% Dataset 'precipitationCal'
% Size: 1800x3600
% MaxSize: 1800x3600
% Datatype: H5T_IEEE_F32LE (single)
% ChunkSize: 1800x145
% Filters: deflate(6)
% FillValue: 0.000000
% Attributes:
% 'DimensionNames': 'lon,lat'
% 'Units': 'mm/hr'
% 'units': 'mm/hr'
% 'coordinates': 'lon lat'
% '_FillValue': -9999.900391
% 'CodeMissingValue': '-9999.9'
% 'DIMENSION_LIST': H5T_VLEN
% Dataset 'precipitationUncal'
% Size: 1800x3600
% MaxSize: 1800x3600
% Datatype: H5T_IEEE_F32LE (single)
% ChunkSize: 1800x145
% Filters: deflate(6)
% FillValue: 0.000000
% Attributes:
% 'DimensionNames': 'lon,lat'
% 'Units': 'mm/hr'
% 'units': 'mm/hr'
% 'coordinates': 'lon lat'
% '_FillValue': -9999.900391
% 'CodeMissingValue': '-9999.9'
% 'DIMENSION_LIST': H5T_VLEN
- 代码解释以及读取数据时的注意事项
% Detailed explanation goes here
% read GPM IMERG preciptation datasets (cal & uncal)
% Details:
% var_name: precipitationCal & preciptationUncal
% dims: 1800 * 3600
% original extent: 180W~180E 90S~90N. -->注意数据要按行翻转、要转换成0—360N
% origin: southwest. cell size:0.1
% First point left corner: 89.9S,179.9W
% Second point left corner: 89.9S,179.8W
% Last point left corner: 89.9N,179.9E
% function: h5disp(fname);
% h5info('...');
% h5read(fname, 'dataset's path');
- 读取代码（Matlab）
这里将读取单个gpm文件的代码封装为一个function, 该函数读取了gpm文件中的未校正（precipitationUncal）和有站点校正（precipitationCal）的数据集(读者可自行修改所需要读取的数据集)。
在文章最后作者会给出一个matlab脚本，来调用这个read_IMERG函数，从而进行gpm数据的批量处理。
function [preUncal, preCal] = read_IMERG( fileFullName )
% 读取IMERG降水数据 (lat, lon) = (1800, 3600)
preUncal = h5read(fileFullName, '/Grid/precipitationUncal'); % 未校正降雨数据
preCal = h5read(fileFullName, '/Grid/precipitationCal'); % 有校正的降雨数据
preUncal = flipdim(preUncal, 1); % 转换成 90N~90S <--> filpud(preUncal);
preCal = flipdim(preCal, 1);
preUncal = [preUncal(:,1801:end), preUncal(:,1:1800)]; % 转换成0—360N
preCal = [preCal(:,1801:end), preCal(:,1:1800)];
preUncal = preUncal ./ 2; % 单位转换成 mm/0.5 hour
preCal = preCal ./ 2;
preUncal(preUncal<0) = -999; % 无效值替换
preCal(preCal<0) = -999;
end
GPM-GSMaP数据
- 数据介绍
-
GSMaP依据不同的传感器输入源和反演算法提供了三套主流类型的卫星降水数据，分别为实时卫星降水数据GSMaP_NRT,标准研究型卫星降水数据GSMaP_MVK,地面校正卫星降水数据GSMaP_Gauge。
GSMaP数据反演流程
GSMaP数据格式
-
数据下载(ftp下载)
数据格式和下载地址
- 数据读取(matlab)
% % gsmap:
% 60N-->60S, 0-->360E;
% 0.1*0.1degree;
% 1 hourly;
% 4-byte-float;
function [rain] = read_GsMap( fullPathName )
% 读取gsmap降水数据 (lon, lat) = (3600,1200)
fid = fopen(fullPathName, 'r');
rain = fread(fid, [3600,1200], 'float');
rain(rain < 0) = -999;
fclose(fid);
end
%---------------------Detailed explanation goes here-----------------------
% * GrADS control file for GSMaP_MVK Hourly Gauge-calibrated Rain (ver.7).
% *
% * Negative value indicates missing due to following reason.
% * -4: missing due to sea ice in microwave retrieval
% * -8: missing due to low temperature in microwave retrieval
% * -99: missing due to no observation by IR and/or microwave
% *
% DSET ^gsmap_gauge.%y4%m2%d2.%h200.v7.0000.0.dat
% TITLE GSMaP_GAUGE 0.1deg Hourly (ver.7)
% OPTIONS YREV LITTLE_ENDIAN TEMPLATE
% UNDEF -99.0
% XDEF 3600 LINEAR 0.05 0.1
% YDEF 1200 LINEAR -59.95 0.1
% ZDEF 1 LEVELS 1013
% TDEF 87600 LINEAR 00Z1jan2014 1hr
% VARS 1
% precip 0 99 hourly averaged rain rate [mm/hr]
% ENDVARS
总结
matlab脚本批量处理GPM卫星降水数据
这里以gpm-imerg数据为例，通过调用上面给出的read_IMERG(filename)函数实现数据的批量读取。从而计算全球区域上某个时间段内的平均降水。gsmap批量处理的方式类似！
读者只需修改数据文件的路径即可！
addpath('E:\work\matlab处理卫星降水\code\fun_lib\', '-end');
clc;
clear;
tic; %计时
% 定义头文件
header = ['ncols 3600\r\n'...
'nrows 1800\r\n'...
'xllcorner 0.00\r\n'...
'yllcorner -90.00\r\n'...
'cellsize 0.1\r\n'...
'NODATA_value -999.00\r\n'];
folder = uigetdir('E:\work\', 'Select a Dir to Open');
precipUncal = zeros(1800,3600);
precipCal = zeros(1800,3600);
valid = zeros(1800,3600); %有效降雨记录
Files = dir(fullfile(folder, '\*.RT-H5'));
if ~isempty(Files)
for i=1:length(Files) %循环单个降雨数据文件
fullpath = fullfile( folder, Files(i).name );
disp(Files(i).name);
[preUncal, preCal] = read_IMERG(fullpath); % 调用自定义函数
%累计降雨
valid(preUncal>=0) = valid(preUncal>=0) + 1;
preUncal(preUncal<0) = 0;
precipUncal = precipUncal + preUncal;
precipCal = precipCal + preCal;
end
end
%save
dailyPre = precipUncal ./ valid .* 2 .* 24; %单位转换 mm/daily
dailyPre(~isfinite(dailyPre)) = -999;
ofid = fopen('E:\work\GPM\data\testIMERG_early.txt', 'w');
fprintf(ofid, header);
fprintf(ofid, strcat(repmat('%8.3f ',[1, 3600]),'\r\n'), dailyPre');
fclose(ofid);
toc;
disp('.............数据处理完成............Good Luck!!!');
Tips:如有任何疑问和不足，欢迎关注并对文章进行评论，作者将进行细心的解答和补充！！！
