常用水文气象数据读取及其可视化（二进制、HDF5、NetCDF）

2020-05-18 本文已影响0人 zengsk

前言

地学、水文、气象领域的自然科学数据通常以netcdf、hdf、二进制等方式存储，比如温度、降水、蒸发数据等；学会这些数据格式的读取和可视化是进行地学统计分析计算的关键，python提供了解析nc、hdf等数据格式的库，作者整理了如何利用python进行.nc、.hdf、*.dat格式数据的读取并进行栅格转换，欢迎转载和关注！

arcpy读取二进制文件---以GSMaP降水数据集为例

温度、降水、蒸发等科学数据为了达到压缩文件大小，易于存储和维护等目的，通常在下载后往往以二进制.bin、.dat等方式呈现。因此在进行数据分析和统计计算时需要读取二进制文件进行数据还原。

注意：读取二进制文件要特别要注意到二进制文件数据存储的端序问题（字节序），小端序需要进行端序转换！

本地文件夹下存在2018年6月份前10天的日尺度降水数据集，下面利用arcoy脚本将该数据批量转换为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小时尺度降水数据，输出为tiff
    2.运行环境需要安装python2 需要arcpy模块
    3.使用arcgis自带的python环境(有arcpy模块)
    4.运行结果可以直接用arcgis打开
'''

# 导入module
import os
import glob
import numpy as np
import arcpy
import warnings

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

# 原始降水数据文件夹,可根据自己本地情况修改
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' # 输出文件名(可根据实际情况改)
    # arcpy.NumPyArrayToRaster()不清楚输入参数可以查看arcpy的官方文档
    # 矩阵转为栅格
    raster = arcpy.NumPyArrayToRaster(ds, arcpy.Point(0, -60.0),
                    x_cell_size=0.1, y_cell_size=0.1, value_to_nodata=-999.00)
    # 添加地理坐标系 GCS_WGS_1984
    spatialRef = arcpy.SpatialReference(4326)
    arcpy.DefineProjection_management(raster, spatialRef)
    raster.save(TiffName)

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

处理结果

image

arcmap中展示

image

NetCDF（*.nc）文件转GeoTIFF---以GLDAS路面蒸发数据为例

在地学领域，nc格式的文件可谓随处可见，netCDF文件可以存储多维数字矩阵，同时又封装了自描述信息（例如经纬度、高度层、时间戳、单位等）。nc文件的IO接口也很普及，Python、NCL、Matlab等气象上常用的软件都可以对其进行读写操作。用Python对nc文件的读写操作是使用Python进行气象数据分析最基础的部分之一。

基于Python处理nc数据必要的库是netCDF4，同时如果需要将nc文件转换为tiff文件，还需要osgeo库中的gdal子库和osr子库。**
下面以GLDAS全球陆地蒸发数据为例，展示如何利用python读取nc文件并将数据转换为geotiff文件。

源码

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

"""
Details:   读取GLDAS_NOAH025的蒸发数据（netcdf格式）并转为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文件的变量名称， 默认 Evap_tavg
    :param save_name:  输出tiff文件的保存路径
    :return:
    """
    print(nc_url)
    nc_fp = nc.Dataset(nc_url)
    # 查看属性
    print(nc_fp)
    # # 查看变量
    print(nc_fp.variables)
    # 获取var数据集
    evap = nc_fp.variables[var][:]  # 单位转换 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)


# 数据矩阵转tiff文件
def array2raster(array, xsize, ysize, rasterOrigin, newTiffName):
    """
     array: 计算后的栅格数据
     xsize: x方向像元大小
     ysize: y方向像元大小
     rasterOrigin: 原始栅格数据Extent
     newRasterfn: 输出tif路径
    """
    cols = array.shape[1]  # 矩阵列数
    rows = array.shape[0]  # 矩阵行数
    originX = rasterOrigin[0]  # 起始像元经度
    originY = rasterOrigin[1]  # 起始像元纬度
    driver = gdal.GetDriverByName('GTiff')
    outRaster = driver.Create(newTiffName, cols, rows, 1, gdal.GDT_Float32)
    # 设置仿射矩阵参数（左上角X坐标; X方向分辨率; 旋转角度，如果图像北方朝上，该值为0; 左上角Y坐标; 旋转角度; y方向分辨率）
    outRaster.SetGeoTransform((originX, xsize, 0, originY, 0, ysize))
    # 获取数据集第一个波段，是从1开始，不是从0开始
    outband = outRaster.GetRasterBand(1)
    outband.WriteArray(array)
    spatialRef = osr.SpatialReference()
    # 代码4326表示WGS84坐标
    spatialRef.ImportFromEPSG(4326)
    outRaster.SetProjection(spatialRef.ExportToWkt())
    outband.FlushCache()
    outRaster.FlushCache()


# 主函数
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 ! ++++++")

结果在arcmap中显示

image

读取hdf文件（*.hdf）并转为tif文件---以MODIS地面温度数据为例

HDF文件是地学研究中常用的数据格式，卫星数据的储存格式通常如此。下面以MODIS 16A3 地表温度数据为例，进行LST数据的读取和TIFF格式转换。

利用python读取hdf文件需要用的pyhdf库，lst数据转换为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  # 转换为摄氏度
    array2raster(lst, 1.0, -1.0, [-180, 90], save_name)


def array2raster(array, xsize, ysize, rasterOrigin, newTiffName):
    """
     array: 计算后的栅格数据
     xsize: x方向像元大小
     ysize: y方向像元大小
     rasterOrigin: 原始栅格数据Extent
     newRasterfn: 输出tif路径
    """
    cols = array.shape[1]  # 矩阵列数
    rows = array.shape[0]  # 矩阵行数
    originX = rasterOrigin[0]  # 起始像元经度
    originY = rasterOrigin[1]  # 起始像元纬度
    driver = gdal.GetDriverByName('GTiff')
    outRaster = driver.Create(newTiffName, cols, rows, 1, gdal.GDT_Float32)
    # 设置仿射矩阵参数（左上角X坐标; X方向分辨率; 旋转角度，如果图像北方朝上，该值为0; 左上角Y坐标; 旋转角度; y方向分辨率）
    outRaster.SetGeoTransform((originX, xsize, 0, originY, 0, ysize))
    # 获取数据集第一个波段，是从1开始，不是从0开始
    outband = outRaster.GetRasterBand(1)
    outband.WriteArray(array)
    spatialRef = osr.SpatialReference()
    # 代码4326表示WGS84坐标
    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 ! ++++++")