常用水文氣象數(shù)據(jù)讀取及其可視化（二進制、HDF5、NetCDF）--以GLDAS、MODIS、GSMaP為例

前言

地學(xué)、水文、氣象領(lǐng)域的自然科學(xué)數(shù)據(jù)通常以netcdf、hdf、二進制等方式存儲，比如溫度、降水、蒸發(fā)數(shù)據(jù)等；學(xué)會這些數(shù)據(jù)格式的讀取和可視化是進行地學(xué)統(tǒng)計分析計算的關(guān)鍵，python提供了解析nc、hdf等數(shù)據(jù)格式的庫，作者整理了如何利用python進行.nc、.hdf、*.dat格式數(shù)據(jù)的讀取并進行柵格轉(zhuǎn)換，歡迎轉(zhuǎn)載和關(guān)注！

arcpy讀取二進制文件---以GSMaP降水?dāng)?shù)據(jù)集為例

溫度、降水、蒸發(fā)等科學(xué)數(shù)據(jù)為了達到壓縮文件大小，易于存儲和維護等目的，通常在下載后往往以二進制.bin、.dat等方式呈現(xiàn)。因此在進行數(shù)據(jù)分析和統(tǒng)計計算時需要讀取二進制文件進行數(shù)據(jù)還原。

注意： 讀取二進制文件要特別要注意到 二進制文件數(shù)據(jù)存儲的端序問題（字節(jié)序），小端序需要進行端序轉(zhuǎn)換！

  • 本地文件夾下存在2018年6月份前10天的日尺度降水?dāng)?shù)據(jù)集，下面利用arcoy腳本將該數(shù)據(jù)批量轉(zhuǎn)換為tiff柵格影像，并在arcmap中展示。
image

代碼展示

#!/usr/bin/env python
# -*- coding:utf-8 -*-
# Name   : gsmap2TIFF.py
# Author : zengsk in NanJing
# Created: 2019/8/24 23:26

'''
說明：1.該腳本是讀取gsmap小時尺度降水?dāng)?shù)據(jù)，輸出為tiff
    2.運行環(huán)境需要安裝python2 需要arcpy模塊
    3.使用arcgis自帶的python環(huán)境(有arcpy模塊)
    4.運行結(jié)果可以直接用arcgis打開
'''

# 導(dǎo)入module
import os
import glob
import numpy as np
import arcpy
import warnings

warnings.simplefilter("ignore") # 忽略警告

# 原始降水?dāng)?shù)據(jù)文件夾,可根據(jù)自己本地情況修改
sPath = r'H:\test\gsmap'
oDir = r"H:\test\res"

for fileName in glob.glob(sPath+'\*.dat'):
    print("Processing... {0}".format(fileName))
    ds = np.fromfile(fileName, dtype=np.float32)
    ds = np.resize(ds,[1200,3600])
    ds[ds < 0] = -999.00
    ds[np.isnan(ds)] = -999.00 # NODATA_value
    
    # 輸出為TIFF(注意：要用到arcpy模塊)
    if not os.path.exists(oDir):
        os.makedirs(oDir)
    TiffName = oDir + os.sep + os.path.basename(fileName)[0:-4] + '.tif' # 輸出文件名(可根據(jù)實際情況改)
    # arcpy.NumPyArrayToRaster()不清楚輸入?yún)?shù)可以查看arcpy的官方文檔
    # 矩陣轉(zhuǎn)為柵格
    raster = arcpy.NumPyArrayToRaster(ds, arcpy.Point(0, -60.0),
                    x_cell_size=0.1, y_cell_size=0.1, value_to_nodata=-999.00)
    # 添加地理坐標(biāo)系 GCS_WGS_1984
    spatialRef = arcpy.SpatialReference(4326)
    arcpy.DefineProjection_management(raster, spatialRef)
    raster.save(TiffName)

print("\n++++++ Data Processing Successfully Completed ! ++++++")

處理結(jié)果

image

arcmap中展示

image

NetCDF（*.nc）文件轉(zhuǎn)GeoTIFF---以GLDAS路面蒸發(fā)數(shù)據(jù)為例

在地學(xué)領(lǐng)域，nc格式的文件可謂隨處可見，netCDF文件可以存儲多維數(shù)字矩陣，同時又封裝了自描述信息（例如經(jīng)緯度、高度層、時間戳、單位等）。nc文件的IO接口也很普及，Python、NCL、Matlab等氣象上常用的軟件都可以對其進行讀寫操作。用Python對nc文件的讀寫操作是使用Python進行氣象數(shù)據(jù)分析最基礎(chǔ)的部分之一。

  • 基于Python處理nc數(shù)據(jù)必要的庫是netCDF4，同時如果需要將nc文件轉(zhuǎn)換為tiff文件，還需要osgeo庫中的gdal子庫和osr子庫。**

  • 下面以GLDAS全球陸地蒸發(fā)數(shù)據(jù)為例，展示如何利用python讀取nc文件并將數(shù)據(jù)轉(zhuǎn)換為geotiff文件。

源碼

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Name   : gldas2tiff.py
# Author : zengsk in NanJing
# Created: 2020/5/13 17:02

"""
Details:   讀取GLDAS_NOAH025的蒸發(fā)數(shù)據(jù)（netcdf格式）并轉(zhuǎn)為Tiff柵格影像
"""

import netCDF4 as nc
from osgeo import gdal, osr
import numpy as np


# ('Evap_tavg', <class 'netCDF4._netCDF4.Variable'>
# float32 Evap_tavg(time, lat, lon)
#     units: kg m-2 s-1
#     standard_name: water_evaporation_flux
#     long_name: Evapotranspiration
#     cell_methods: time: mean
#     scale_factor: 1.0
#     add_offset: 0.0
#     missing_value: -9999.0
#     _FillValue: -9999.0
#     vmin: -1.8462154e-06
#     vmax: 7.190982e-05
# unlimited dimensions: time
# current shape = (1, 600, 1440)
# filling off
# )


def nc2tiff(nc_url, var='Evap_tavg', save_name='./nc2tiff.tif'):
    """
    :param nc_url:  netcdf文件的路徑
    :param var:   讀取nc文件的變量名稱， 默認(rèn) Evap_tavg
    :param save_name:  輸出tiff文件的保存路徑
    :return:
    """
    print(nc_url)
    nc_fp = nc.Dataset(nc_url)
    # 查看屬性
    print(nc_fp)
    # # 查看變量
    print(nc_fp.variables)
    # 獲取var數(shù)據(jù)集
    evap = nc_fp.variables[var][:]  # 單位轉(zhuǎn)換 kg m-2 s-1
    evap = np.array(evap)
    evap = np.squeeze(evap)
    # missing value 處理
    fillvalue = nc_fp.variables[var].missing_value
    evap[evap == fillvalue] = np.nan
    evap = evap * 3600 * 24 * 30
    array2raster(evap, 0.25, 0.25, [-180, -60], save_name)


# 數(shù)據(jù)矩陣轉(zhuǎn)tiff文件
def array2raster(array, xsize, ysize, rasterOrigin, newTiffName):
    """
     array: 計算后的柵格數(shù)據(jù)
     xsize: x方向像元大小
     ysize: y方向像元大小
     rasterOrigin: 原始柵格數(shù)據(jù)Extent
     newRasterfn: 輸出tif路徑
    """
    cols = array.shape[1]  # 矩陣列數(shù)
    rows = array.shape[0]  # 矩陣行數(shù)
    originX = rasterOrigin[0]  # 起始像元經(jīng)度
    originY = rasterOrigin[1]  # 起始像元緯度
    driver = gdal.GetDriverByName('GTiff')
    outRaster = driver.Create(newTiffName, cols, rows, 1, gdal.GDT_Float32)
    # 設(shè)置仿射矩陣參數(shù)（左上角X坐標(biāo); X方向分辨率; 旋轉(zhuǎn)角度，如果圖像北方朝上，該值為0; 左上角Y坐標(biāo); 旋轉(zhuǎn)角度; y方向分辨率）
    outRaster.SetGeoTransform((originX, xsize, 0, originY, 0, ysize))
    # 獲取數(shù)據(jù)集第一個波段，是從1開始，不是從0開始
    outband = outRaster.GetRasterBand(1)
    outband.WriteArray(array)
    spatialRef = osr.SpatialReference()
    # 代碼4326表示W(wǎng)GS84坐標(biāo)
    spatialRef.ImportFromEPSG(4326)
    outRaster.SetProjection(spatialRef.ExportToWkt())
    outband.FlushCache()
    outRaster.FlushCache()


# 主函數(shù)
if __name__ == '__main__':
    netcdf_url = r'E:\Scripts\Test\GLDAS_NOAH025_M.A200002.021.nc4'
    nc2tiff(netcdf_url, "Evap_tavg", save_name='E:\Scripts\Test\evap2tiff.tif')

    print("\n++++++ Data Processing Successfully Completed ! ++++++")

結(jié)果在arcmap中顯示

image

讀取hdf文件（*.hdf）并轉(zhuǎn)為tif文件---以MODIS地面溫度數(shù)據(jù)為例

HDF文件是地學(xué)研究中常用的數(shù)據(jù)格式，衛(wèi)星數(shù)據(jù)的儲存格式通常如此。下面以MODIS 16A3 地表溫度數(shù)據(jù)為例，進行LST數(shù)據(jù)的讀取和TIFF格式轉(zhuǎn)換。

  • 利用python讀取hdf文件需要用的pyhdf庫，lst數(shù)據(jù)轉(zhuǎn)換為geotiff文件則采用GDAL庫。

源碼

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Name   : modis2tiff.py
# Author : zengsk in NanJing
# Created: 2020/5/17 20:28

"""
Details:
"""

from pyhdf import SD
from osgeo import gdal, osr
import numpy as np


def hdf2tiff(hdf_url, variable, save_name="./hdf2tiff.tiff"):
    fp = SD.SD(hdf_url)
    for idx, sds in enumerate(fp.datasets().keys()):
        print(idx, sds)
    # print(fp.select(variable).attributes())
    lst = fp.select(variable).get()
    lst[lst == -9999.0] = np.nan
    lst -= 273.15  # 轉(zhuǎn)換為攝氏度
    array2raster(lst, 1.0, -1.0, [-180, 90], save_name)


def array2raster(array, xsize, ysize, rasterOrigin, newTiffName):
    """
     array: 計算后的柵格數(shù)據(jù)
     xsize: x方向像元大小
     ysize: y方向像元大小
     rasterOrigin: 原始柵格數(shù)據(jù)Extent
     newRasterfn: 輸出tif路徑
    """
    cols = array.shape[1]  # 矩陣列數(shù)
    rows = array.shape[0]  # 矩陣行數(shù)
    originX = rasterOrigin[0]  # 起始像元經(jīng)度
    originY = rasterOrigin[1]  # 起始像元緯度
    driver = gdal.GetDriverByName('GTiff')
    outRaster = driver.Create(newTiffName, cols, rows, 1, gdal.GDT_Float32)
    # 設(shè)置仿射矩陣參數(shù)（左上角X坐標(biāo); X方向分辨率; 旋轉(zhuǎn)角度，如果圖像北方朝上，該值為0; 左上角Y坐標(biāo); 旋轉(zhuǎn)角度; y方向分辨率）
    outRaster.SetGeoTransform((originX, xsize, 0, originY, 0, ysize))
    # 獲取數(shù)據(jù)集第一個波段，是從1開始，不是從0開始
    outband = outRaster.GetRasterBand(1)
    outband.WriteArray(array)
    spatialRef = osr.SpatialReference()
    # 代碼4326表示W(wǎng)GS84坐標(biāo)
    spatialRef.ImportFromEPSG(4326)
    outRaster.SetProjection(spatialRef.ExportToWkt())
    outband.FlushCache()


if __name__ == '__main__':
    hdf_url = r'E:\Scripts\Test\LST\MOD11CM1_v005_day___lst_2000m03.hdf'
    hdf2tiff(hdf_url, "day_lst", save_name=r'E:\Scripts\Test\LST_TIFF2.tif')
    print("\n++++++ Data Processing Successfully Completed ! ++++++")
image

處理結(jié)果在arcmap中顯示

image
最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容