基于Python gdal包的影像投影及重采样

动机

在做论文的时候，使用到了AMSR-2的数据，其原始数据为.h5格式，需要“FormatConversionTool.exe”软件批量转换为.tif文件。问题来了，转换之后的.tif文件没有没有投影坐标系，且没有设置Nodata值(图1)，黑乎乎的一片。最大值最小值相等的情况下，可能会误导初使用者，认为数据有问题（其实没问题的啦，只是这个.tif文件没有统计信息而已，你只需要数据右键，重新导出，arcgis会自动计算统计值）。

• 图1

地理相关数据，拿到手，首先要做的就是统一坐标信息。我的研究区是古尔班通古特沙漠，UTM_45度带，因此我需要将自己的数据进行投影。AMSR-2的L1R产品，官方说的是重采样间隔约10km，所以嘛我也就将我的数据投影的同时，进行重采样处理，同时将65535设置为Nodata。因为要处理的数据太多了，所以使用实现（ 其实这些功能arcgis都有，通过arcpy都能够实现的，但是arcpy的效率好像稍微的低一点）。

AMSR-2 L1R 数据简介

L1R 数据是轨道被动微波亮温数据，每个扫描轨道为一个.h5文件（如GW1AM2_2018010105 53_231A_L1 SGRTBR_ 2220220.h5），里面包含了许多信息（不同通道，同一通道又分为垂直极化与水平极化，相同通道及极化方式，又分为不同的采样间隔），具体自己下载一个看看就好。图2给出了L1R .h5文件的部分数据集。因为数据较为复杂，建议初接触AMSR-2数据的去下载L3产品，L3产品将每日全球的全部合成在一张图上了（如果研究区不覆盖南北极，不需要下载文件名中带PN与PS产品，只下载EQ就好）。详细信息去看JAXA的G-Portal网站。

• 图2

FormatConversionTool

FormatConversionTool软件，可以将所有的AMSR-2产品，转化为.tif文件，速度也挺快的，下载地址依然是在Jaxa的G-Portal，简单好用。但是我一直希望有一个python版本，直接将.h5文件转化为.tif，这样我的我的所有处理数据的代码就能集成起来了，但是没找到python版本(如果谁找到，记得跟我说下)。FormatConversionTool软件，默认设置，会将所有的AMSR-2 L1R的.h5文件自己转换为.tif文件，如果只需要部分子集，请取消勾选，如图3。转换此文件转换完成，就变成了图1的样子，需要投影了

• 图3

基于 Python—Gdal—Warp的投影转换

因为数据量有点多，我想估计下时间，所以我加了个记录时间的代码；其他代码中标注的挺仔细的，所以我就不解释。

#!/usr/bin/python
# -*- coding: utf-8 -*-

"""
@version: Anaconda
@author: LeYongkang
@contact: 1363989042@qq.com
@software: PyCharm
@file: 01_AMSR2_Prj_Gdal.py
@time: 2020/3/4 0004 下午 7:22
"""


"""
NDVI EVI AMSR2 LST

重投影时，重采样

重投影，投影过程中可以设置以下信息
gdal.Warp(
    xRes,yRes:两个方向上的分辨率；
    srcNodata：原来数据的Nodata值
    dstNodata：输出数据的Nodata值
    dstSRS：输出的投影坐标系，可以读取影像的：
            Raster = gdal.Open(InputImage, gdal.GA_ReadOnly)
            Projection = Raster.GetProjectionRef()
    resampleAlg:重采样方式,算法包括：
            import gdalconst
            gdalconst.GRA_NearestNeighbour：near
            gdalconst.GRA_Bilinear:bilinear
            gdalconst.GRA_Cubic:cubic
            gdalconst.GRA_CubicSpline:cubicspline
            gdalconst.GRA_Lanczos:lanczos
            gdalconst.GRA_Average:average
            gdalconst.GRA_Mode:mode
    )

"""
import numpy as np
import os,datetime,gdal,gdalconst



if __name__ == "__main__":

    # 获取程序开始时间
    start = datetime.datetime.now()
    # Open datasets，此处的tif文件为模板文件，其提供目标投影坐标
    InputImage = r'F:\Parper2\Tif_Mask\Mask_sub10km\AMSR2_TifMask_UTM_45.tif'
    Raster = gdal.Open(InputImage, gdal.GA_ReadOnly)
    Projection = Raster.GetProjectionRef()

    # 需要重投影的tif文件所在文件夹
    inDir = r"I:\2018_Parper2\AMSR2_2018\.tif"
    # 投影后tif文件存放文件夹
    outDir = r"I:\2018_Parper2\AMSR2_2018\.GeoTif"

    # 如果没有输出文件夹，创建文件夹
    if not os.path.exists(outDir):
        os.makedirs(outDir)

    # 获取文件，同时剔除 "00" 的文件（含有00的为月均）
    inList = [name for name in os.listdir(inDir) if name.endswith(".tif") and "00" not in name.split("_")[1]]
    print(inList)

    for file in inList:
        # 输入文件的完整路径
        input_raster = os.path.join(inDir,file)
        # 输出文件的完整路径
        output_raster = os.path.join(outDir,file)
        if not os.path.exists(output_raster):
            OutTile = gdal.Warp(output_raster, input_raster,
                      dstSRS=Projection,
                      # cutlineDSName=r'F:\Parper2\Shp_Mask\RectangleOutline.shp',
                      resampleAlg=gdalconst.GRA_NearestNeighbour,
                      dstNodata=0,
                      srcNodata=65535,
                      xRes=10000,
                      yRes=10000
                      )
            OutTile =None
            print(inList.index(file))

    # 获取程序结束时间
    end = datetime.datetime.now()
    # 获取共使用了多长时间
    SpendTimne = end - start
    # 新建 .txt文件，写入程序运行时间信息
    with open(outDir + os.sep + 'TimeInformation.txt', 'a') as file_handle:
        file_handle.write(str(start)+'\n')
        file_handle.write(str(end)+'\n')
        file_handle.write(str(SpendTimne)+'\n')
        # 有时放在循环里面需要自动转行，不然会覆盖上一条数据
        file_handle.write('\n')

大家感觉Python的地理处理稍微有些复杂，可以通过Arcpy入门的，等爱上代码的时候，再尝试其他包也未尝不可。

结语

这个博客文章，其实很简单，就是GIS的一个小功能，希望大家也能入手代码，这样以后遇见大量数据，也不会被吓到，导师开心，你也开心！！！！后续可能会仔细说下AMSR-2的L1R 数据如何使用Python下载、机器学习、云服务器等，反正就是我已经做过的工作吧！