EASE-Grid投影风云卫星产品地理信息写入方法

韩书新,安英玉,高 昂,于 敏,秦 铁,王志晓

(1.黑龙江省生态气象中心,黑龙江 哈尔滨 150030;2.黑龙江省人民政府人工降雨办公室,黑龙江 哈尔滨 150030;3.国家卫星气象中心,北京 100081;4.黑龙江省气象局机关服务中心,黑龙江 哈尔滨 150030;5.西安理工大学 计算机科学与工程学院,陕西 西安 710048)

0 引 言

随着国内卫星工程的迅速发展及气象卫星相关技术如晴空数据合成[1]、地面测距精度和卫星数据可视化等卫星应用技术能力的提高[2-5],国内气象卫星遥感应用能力得到迅速提升。目前,气象卫星在林火、水体洪涝、植被、降水、积雪等陆表监测方面[6-10]及天气预报、台风等大气环境监测方面[11-13]的应用提供重要的数据产品支撑;同时,作为全球地球观测系统(GEOSS)的重要成员[14],风云卫星积极参与国际多卫星集成业务系统的研发[15-16];中国在气象卫星应用效益评估方法上也进行了深入研究[17]。

国家卫星气象中心在风云卫星遥感数据服务网发布了很多卫星遥感应用产品数据集,得益于国家卫星气象中心科学的反演算法和高效的数据处理能力,这些产品数据集可实现实时、批量生产[18-20]。对于省级卫星遥感应用来说,采用这些国家级产品数据集对提高省级遥感业务能力、标准化程度和工作效率都有着积极的作用。

目前,风云卫星产品数据集在黑龙江省的利用率较低。由于省级卫星遥感能力上的不足,导致在数据集的使用中也存在问题。该文通过基于EASE-Grid投影的风云系列卫星遥感产品数据集在省级应用中发现的问题进行研究并提出解决方法,对提高风云卫星遥感产品的省级应用和提升省级卫星遥感业务能力水平都有着积极的作用。

1 EASE-Grid投影简介

等面积可伸缩地球网格(Equal-Area Scalable Earth Grid,EASE-Grid),是依据美国雪冰数据中心(National Snow and Ice Data Center,NSIDC)的数字地图和网格定义理论设计的一种网格,假定网格数据集被完全定义为地图投影和网格点的覆盖栅格,最初用于美国国家海洋和大气管理局(National Oceanic and Atmospheric Administration,NOAA)和美国国家航空航天局(National Aeronautics and Space Administration,NASA)的开拓者计划的微波成像(SSM/I)项目制作数据产品开发出来的[21]。

EASE-Grid包含3种投影,采用其投影的数据可被较好的扩展和应用,可作为通用工具用于处理全球尺度的网格数据。作为全球尺度网格数据(包括处理遥感数据)的通用投影格式,EASE-Grid投影格式是一种等面积投影,可使数据表示为具有多种网格分辨率的数字阵列(12.5 km和25 km等),依据适用范围的不同分为3种子投影：全球圆柱等面积投影(Global)、北半球方位角等面积投影(Northern Hemisphere)、南半球方位角等面积投影(Southern Hemisphere)。都是基于“1984全球大地坐标系统(World Geodetic System 1984,WGS84)”的椭球体,但投影方式有所差异,其中,EASE-Grid的Global属于“等面积圆柱投影(Cylindrical Equal Area)”,属于一种切投影,标准纬线为赤道;“EASE-Grid Northern Hemisphere”和“EASE-Grid Southern Hemisphere”采用“兰伯特等面积方位投影(Lambert Azimuthal Equal Area)”,切点分别为南、北极点[22-23]。

目前EASE-Grid有两个版本,即EASE-Grid和EASE-Grid2.0,中国风云气象卫星遥感产品设计中主要使用的是EASE-Grid版,该文主要对EASE-Grid的数据格式进行分析处理。EASE-Grid投影区域如图1所示。

图1 EASE-Grid投影区域

EASE-Grid地图投影包括：南、北半球方位角投影,全球圆柱投影和温带圆柱投影(不含灰色阴影区域)。其中,温带圆柱投影属于EASE-Grid2.0投影,覆盖的纬度范围为(67°S,67°N),在其覆盖范围内与全球圆柱投影相同,全球圆柱投影纬度范围为(84°S,84°N)。25 km分辨率的EASE-Grid参数如表1所示。

表1 25 km分辨率的3种EASE-Grid投影参数

EASE-Grid的3种投影定义[22]如下：

南半球的方位角等面积投影定义如公式1～4所示：

r=2×R/C×sin(lon)×cos(π/4-lat/2)+r0

(1)

s=-2×R/C×cos(lon)×cos(π/4-lat/2)+s0

(2)

h=sin(π/4-lat/2)

(3)

k=csc(π/4-lat/2)

(4)

北半球方位角等面积投影定义如公式5～8所示：

r=2×R/C×sin(lon)×sin(π/4-lat/2)+r0

(5)

s=2×R/C×cos(lon)×sin(π/4-lat/2)+s0

(6)

h=cos(π/4-lat/2)

(7)

k=sec(π/4-lat/2)

(8)

全球圆柱等面积投影定义如公式9～12所示：

r=r0+R/C×lon×cos(30°)

(9)

s=s0-R/C×sin(lat)/cos(30°)

(10)

h=cos(lat)/cos(30°)

(11)

k=cos(30°)/cos(π)

(12)

其中,r为列号,s为行号,h为沿经线方向的比例,k为沿纬线方向的比,lon为经度(弧度),lat为纬度(弧度),R为地球半径(km),C为像元大小(km),r0为地图原点列号,s0为地图原点行号。

事实上，我所付出的辛劳并没有白费，这些法律知识的编译，让旅居异国的朋友们在生活上获得了指标。此后每一年，我都收到很多华人朋友的感谢信函，他们欣喜，我也快慰。助人为快乐之本，谁说不是呢？

2 FY3D雪水当量产品介绍

国家卫星气象中心的FY3D雪水当量(SWE)产品数据集是基于微波成像仪(MWRI)数据开发的,数据类型为HDF5格式,投影方式为南、北半球方位角投影(Southern Hemisphere,Lambert Azimuthal &Northern Hemisphere,Lambert Azimuthal),分辨率为25 km,数据集包括雪深(SD)雪水当量(SWE)全球区域(南、北半球)升/降轨产品[24-26],数据集描述如表2所示。

用ArcMap导入FY3D的雪水当量日产品,输出图像如图2所示。可以看出,FY3D雪水当量产品的地理坐标不是常规的WGS84坐标系,而是EASE-Grid的北半球方位角投影,并且其数据本身没有附带经纬度等地理信息,无法与基于WGS84坐标系的全球区域矢量(图中左上角的全球矢量数据)进行匹配。

图2 基于EASE-Grid投影的FY3D雪水当量日产品(升轨)

为了便于后续科研和业务工作中使用GIS类软件对数据进行分析处理,需要将其转换为带经纬度等地理信息的GeoTiff格式的数据。

3 数据投影转换

3.1 GeoTiff文件的地理信息目录结构

GeoTiff是一种Tiff6.0文件,它继承了Tiff6.0文件规范中的相应部分[27],所有的GeoTiff特有的信息都编码在Tiff的一些预留标签(Tag)中,它没有自己的图像文件目录、二进制结构等一些对Tiff来说用来描述GeoTiff不可见的投影参数及类型信息[28]。独立的信息标签会导致标签用量占用过高,进而消耗Tiff有限的标签资源。

为了解决上述问题,GeoTiff文件采用键(Keys)来存储这些信息,这些键在功能上相当于标签,但它处在Tiff的更上一层,与格式化的标签值一起共存,用来支持Tiff文件中的图像数据。这些键也称为GeoKeys,所有键都由‘GeoKeyDirectoryTag’标签来索引,该标签就相当于表示数据地理信息键的一个目录。

3.2 地理信息目录写入方法

国家卫星气象中心风云卫星遥感数据服务网上发布了大量的卫星遥感应用产品数据集。很多产品数据集投影方式采用的都是EASE-Grid的投影方式。下面以FY3D雪水当量日产品北半球投影数据(SWE_Northern_Daily)为例,对原始EASE-Grid投影数据产品的地理信息写入方法进行阐述,并将HDF格式转为GeoTiff格式输出。

通过EASE-Grid参数创建地理坐标系栅格参考对象(GeographicCellsReference)。创建地理坐标系栅格参考对象的主要参数包括25 km分辨率下的x,y坐标系下的范围大小,取值见表1,均为[-9 036 842.76,9 036 842.76],栅格数据矩阵大小为721×721。通过maprefcells命令建立地理坐标系栅格参考对象R,如表3所示。对于北半球地区,需要单独将变量‘ColumnsStartFrom’和‘RowsStartFrom’的值分别设置为‘north’和‘west’。

表3 地理坐标系栅格参考对象R属性

3.2.2 构建地理信息目录

EASE-Grid投影数据所对应的地理信息目录[29](GeoKeyDirectoryTag)的键值共19个,这些键值描述的是EASE-Grid投影坐标系的固有属性并与分辨率无关,故在设置12.5 km分辨率时可采用与25 km相同的标签设置,地理信息目录结构体的键值内容如表4所示。

表4 GeoKeyDirectoryTag主要键值属性

3.2.3 地理信息目录写入

通过geotiffwrite函数将构建好空间参考系和地理信息目录与数据矩阵一起输出为后缀为tif的GeoTiff文件。写入格式为：‘geotiffwrite(‘filename.tif’,Data_Matrix,R,‘GeoKeyDirectoryTag’,info_Northern)’。其中,filename.tif为所要输出的文件名,Data_Matrix为读入的EASE-Grid投影下的雪水当量产品数据矩阵,其矩阵大小应满足721×721大小,R为构建的空间参考系,info_Northern为构建的GeoKeyDirectoryTag。

经过地理信息目录写入的数据以tif文件格式输出,如图3所示,用ArcMap软件导入后,数据文件已带有EASE-Grid北半球投影的经纬度信息,比较图2可以看出矢量地图可与雪水当量数据进行准确的匹配。

图3 NSIDC_EASE_Grid_North投影坐标系下的雪水当量日产品(升轨)

由于数据写入了地理信息,因此可以方便导入GIS类软件进行遥感数据处理,图4为GCS_WGS_1984地理坐标系下的图像数据。

图4 GCS_WGS_1984地理坐标系下的雪水当量日产品(升轨)

通过比较数据写入前后的信息属性可以看出,写入地理信息后的数据已包含EASE-Grid北半球投影的空间参考,如图5所示。

(a)写入前(HDF5格式)

3.3 数据的投影转换

数据写入地理信息后,虽然可匹配GIS数据框属性中设置的坐标系统(WGS1984),但写入地理信息的数据仍为等面积投影,在某些情况下可能涉及到与其它数据像元的经纬度匹配的问题,所以在进行多源或综合数据分析计算时,需要根据需求对数据进行投影转换,如图6所示。

图6 EASE-Grid转WGS1984投影

将EASE-Grid投影下的数据转换到WGS1984投影坐标系统,如图7所示。从图7(a)可以看出,FY3D的雪水当量产品数据在写入地理信息后仍为基于EASE-Grid的一种方位角等面积投影,在WGS1984坐标系统下像元会发生变形;数据经过临近差值法处理后如图7(b)所示,处理后的数据便于匹配其它数据进行后续的分析计算。图7所示的仅是基于EASE-Grid坐标系统到WGS1984坐标系统的投影转换。工作中可根据实际情况需求,在保证数据精度和可用性的情况下,采用不同的差值方法转换为所需的投影方式。

(a)EASE-Grid(栅格范围不一致)

3.4 方法适用性

以上所述的是以FY3D雪水当量产品为例,阐述了EASE-Grid投影数据产品的地理信息写入方法。采用该方法可以解决其它产品的类似问题。图8所示的是经过文中方法处理后的其它数据集产品。写入地理信息后的数据集产品可与矢量数据进行较好的匹配,满足开展科研与业务工作数据需求。

(a)全球陆表温度(LST)

4 结束语

针对基于EASE-Grid投影的国家级卫星遥感产品在省级应用中发现的数据投影问题,以FY3D的雪水当量(SWE)日产品数据为例,通过建立并写入地理信息目录对数据产品进行地理信息写入,解决了卫星遥感产品省级应用中出现的产品数据集与矢量文件或辅助数据(如地形数据等)在分析处理中的地理信息匹配问题,总结如下：

(1)数据处理过程是首先建立与数据相适应的空间参考信息(投影方式、矩阵大小、分辨率等),其次根据空间参考信息构建地理信息目录,最后将已经构建好的空间参考和地理信息目录与原始HDF5格式的数据矩阵一起写入GeoTiff文件类型并输出。

(2)输出的具有地理信息的产品数据仍是EASE-Grid等面积投影,在其它坐标系下图像会发生变形,为满足数据分析处理中的地理信息匹配需要,需根据实际需求对数据进行投影转换操作,再进行后续的数据处理。

(3)采用文中方法也可以用于处理其它产品的类似问题。目前已在FY3D全球陆表温度数据(LST)、海冰密集度(SIC)等数据集产品上进行了适应性验证,经过处理后的数据集产品适用性较好。

针对基于EASE-Grid投影的风云极轨气象卫星遥感反演产品进行地理信息写入,通过投影转换使其适应科研业务应用,对提升省级卫星遥感应用能力有一定的促进作用。该方法虽然在多种反演产品中做了适用性试验,但在具体业务化过程中仍可能出现各种问题,将在今后研究工作中持续改进。