ICEsat数据介绍、下载、读取、GMT筛选

ICESat（Ice, Cloud,and land Elevation Satellite）即冰、云和陆地高程卫星，其主要有效载荷为地球科学激光测高系统（GLAS），是全球首颗对地观测激光测高卫星，于2003年1月13日成功发射。

该卫星由美国德克萨斯州立大学牵头，联合美国NASA和其他工业合作伙伴共同研制，该星主要科学任务包括：1）测量地球两极的冰盖总量，研究冰盖总量对海平面变化的影响；2）测量云和气溶胶的分布和垂直结构；3）测量海冰变化、陆地地形、植被高度以及冰雪面的特征。

但由于激光器故障以及其他多种原因，该星于2010年2月结束科学任务、停止采集数据。在轨运行七年，ICESat不仅为极地变化研究提供了支持，同时为全球高程变化探测提供了近2亿个激光点数据。

NASA于2003年启动了地球观测系统EOS（Earth Observing System）中的世界上首颗激光测高仪试验卫星，第一代冰、云和陆地高程卫星ICESat（Ice，Cloud and IandElevation Satellite），并于2003年1月13日将其搭载在一架波音德尔塔2型运载火箭上在美国加利福尼亚州范登堡空军基地成功发射，主要任务是监测南极洲和格陵兰冰盖的高程变化，从而估算冰盖表面物质平衡和气候的季节及年际变化率。ICESat作为一颗近地激光测高卫星，相比星载雷达高度计提供了更加精确的观测数据；长达7年左右的观测周期较其他空载、机载激光高度计，开展了更长时间序列的测量任务。在极地冰雪环境遥感应用中，ICESat展现出不可替代的作用，突出地表现在冰盖高程测量、海冰厚度测量方面。此外，在植被及生物量测量方面也充分展现了尺度和精度上的优势。

在气候变化背景下，ICESat通过测量极地冰面高程的变化估算冰盖物质平衡和时空变化建立南北极冰盖对海平面变化贡献及其与气候状况的关系模型。该系统被用于测量冰原地形及其相应的温度变化，同时也可用于监测云层和大气的特性。掌握对其时空变化的驱动因素是建立精准的预测模型的关键，通过模拟能够可靠地估计出冰盖对海平面的贡献值。

ICESat卫星上搭载的地球科学激光测高系统GLAS（Geoscience Laser Altimeter System）是第一个用于连续观测地球的激光系统。GLAS系统由供电电源、参考望远镜、主控电路箱、监控板、观星摄像机、雷达监测和电路组件、热管散射系统、测高仪监测组件和三个激光器、激光光束调节机构等共同组成。

该系统共计携带了三台激光器，其中第一台激光器由于使用含铟焊料过多使金属导体氧化成不导电的金铟化合物，造成二极管阵列内部导体被腐蚀，在仅工作了37天后就停止传送数据；第二台激光器也在开始不久后由于激光器的倍频器出现影像暗化，气态碳氢化合物从激光器的黏合处泄漏并与532nm光子产生反应，导致能量快速消耗而停止采集数据，因此NASA改变第三台激光器的运行模式，从不间断测量转为每年进行91天精确重复轨道周期进行测量。新的周期模式分为三个分周期，每个分周期约为33天，调整后激光器表现出色，5年半内一共采集了15个33天数据，采集激光点近2亿个，直至2009年10月11号GLAS停止采集高程数据。【参考：知乎作者：满天星】
详细介绍：ICESat数据介绍及下载

Matlab读取ICEsat GLA14文件程序：

clc;
clear;
fileFolder=fullfile('E:\sunmingzhi\ICEsat_GLAH14\');
dirOutput=dir(fullfile(fileFolder,'*.H5'));
%fileNames={dirOutput.name};

for i=1:size(dirOutput,1)
    dirOutput(i).fullname=strcat(fileFolder,dirOutput(i).name);
    % set the netCDF filename:
    nc_file=dirOutput(i).fullname;
    
    % Display netCDF file.
    ncdisp(nc_file);
    
    % Read data from netCDF file.时间 高度 经度 纬度 共四个参数
    time = ncread(nc_file,'/Data_40HZ/Time/d_UTCTime_40');
    d_elev = ncread(nc_file,'/Data_40HZ/Elevation_Surfaces/d_elev');
    lat = ncread(nc_file,'/Data_40HZ/Geolocation/d_lat');
    lon = ncread(nc_file,'/Data_40HZ/Geolocation/d_lon');
    
    % Load the needed data into an Array.
    nc_data=[lon,lat,time,d_elev];
    
    % find the NAN in the Array and delete the rows.
    [row,col]=find(isnan(nc_data));
    nc_data(row,:)=[];
    
    % Save the data as ascii file.
    dlmwrite(strcat(dirOutput(i).fullname(1:end-3),'.txt'), nc_data, 'delimiter', '\t', 'precision', '%20.6f');
    
end

数据提取后使用GMT筛选研究区域内数据。

数据格式： 经度 纬度 时间 高度
226.613294        53.153328    229650232.520000       740.870000
226.612901        53.154872    229650232.545000       751.692000
226.612115        53.157969    229650232.595000       675.836000
226.611724        53.159512    229650232.620000       760.045000
226.611332        53.161060    229650232.645000       778.305000
226.610940        53.162607    229650232.670000       815.627000
226.610546        53.164154    229650232.695000       834.013000
226.610151        53.165701    229650232.720000       831.373000
226.609755        53.167240    229650232.745000       856.580000

GMT数据筛选
gmtselect： 筛选符合某个特定准则的数据
该命令会从输入文件中读取前两列作为经度和纬度，并判断该点是否符合一定的空间
准则，以筛选出符合条件的记录。输入文件中仅前两列会被使用。
七个空间准则包括：

1. 在矩形区域内（ -R 和 -J ）
2. 与点文件中的每个点的距离在一定范围之内
3. 与线文件中的每条线的距离在一定范围之内
4. 在多边形文件所指定的多边形内
5. 在某个地理区域内（需要海岸线数据）
6. Z 值在某个范围内
7. 该点所在的网格单元内具有有效值（即非零和非 NaN 的值）
   筛选湖泊内的点

gmtselect2103_002_0393_0_01_0001.txt-Rlonmin/lonmix/latmin/latmax-Dh-Nk/s>2103_002_0393_0_01_0001_select.txt

筛选与 pts.txt 中所有点的距离在 300 km 以内，与 lines.txt 中线段的距离在 100km 以外的点:

gmt select lonlatfile -fg -Cpts.txt+d300k -Llines.txt+d100k -Il > subset

此处需要使用 -fg 以告知程序正在处理地理数据。
筛选某个区域内所有不在陆地上的点:

gmt select data.txt -R120/121/22/24 -Dh -Nk/s > subset

筛选 quakes.txt 中所有位于多边形区域内的点:

gmt select quakes.txt -Flonlatpath.txt -fg > subset

stations.txt 中的点投影后与 origin.txt 的距离在 5 cm 之内的点:

gmt select stations.txt -Corigin.txt+d5 -R20/50/-10/20 -JM20c \PROJ_LENGTH_PROJ_LENGTH_UNIT=cm > subset

未完待续…

致谢：
NSIDC（美国冰雪数据中心）提供的数据
GMT中文社区
小李同学Matlab帮助
泽哥Python帮助

E-mail: [e-mail]