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投影的風云極軌氣象衛星遙感反演產品進行地理信息寫入,通過投影轉換使其適應科研業務應用,對提升省級衛星遙感應用能力有一定的促進作用。該方法雖然在多種反演產品中做了適用性試驗,但在具體業務化過程中仍可能出現各種問題,將在今后研究工作中持續改進。