MODIS AQUA 数据 - 使用 python GDAL 堆叠/马赛克数据

Question

我知道如何使用 gdal 和 python 访问和绘制子数据集。但是，我想知道是否有办法使用 HDF4 文件中包含的 GEO 数据，这样我就可以多年查看同一区域。

如果可能，是否可以从数据中切出一个区域以及如何切出？

更新：

更具体地说：我绘制了 MODIS 数据，正如您在下方看到的河流向下移动（矩形结构左上角）。所以整整一年的时间，我观察的位置都不是同一个。

子数据集中有一个名为 Geolocation Fields 的目录，其中包含 Long 和 Alt 目录。那么是否可以访问这些信息或将其覆盖在数据上以切出特定区域？

例如，如果我们看一下下面的 NASA 图片，是否可以将其剪切到 10-15 alt 之间。和 -5 到 0 长。

您可以通过复制以下网址来下载示例文件：

https://ladsweb.modaps.eosdis.nasa.gov/archive/allData/6/MYD021KM/2009/034/MYD021KM.A2009034.1345.006.2012058160107.hdf

更新：

我跑了

x0, dx, dxdy, y0, dydx, dy = hdf_file.GetGeoTransform()

这给了我以下输出：

x0:  0.0
dx:  1.0
dxdy:  0.0
y0:  0.0
dydx:  0.0
dy:  1.0

还有

gdal.Warp(workdir2+"/output.tif",workdir1+"/MYD021KM.A2009002.1345.006.2012058153105.hdf")

这给了我以下错误：

ERROR 1: Input file /Volumes/Transcend/Master_Thesis/Data/AQUA_002_1345/MYD021KM.A2009002.1345.006.2012058153105.hdf has no raster bands.

**更新 2：**

这是我如何打开和读取 hdf 文件的代码：

all_files 是一个包含文件名的列表，例如：

MYD021KM.A2008002.1345.006.2012066153213.hdf
MYD021KM.A2008018.1345.006.2012066183305.hdf
MYD021KM.A2008034.1345.006.2012067035823.hdf
MYD021KM.A2008050.1345.006.2012067084421.hdf
etc .....


for fe in all_files:
    print "\nopening file: ", fe
    try:
        hdf_file = gdal.Open(workdir1 + "/" + fe)
        print "getting subdatasets..."
        subDatasets = hdf_file.GetSubDatasets()
        Emissiv_Bands = gdal.Open(subDatasets[2][0])
        print "getting bands..."
        Bands = Emissiv_Bands.ReadAsArray()
        print "unit conversion ... "

        get_name_tag = re.findall(".A(\d{7}).", all_files[i])[0]
        print "name tag of current file: ", get_name_tag

        # Code for 1 Band:
        L_B_1 = radiance_scales[specific_band] * (Bands[specific_band] - radiance_offsets[specific_band])  # Source: MODIS Level 1B Product User's Guide Page 36 MOD_PR02 V6.1.12 (TERRA)/V6.1.15 (AQUA)
        data_1_band['%s' % get_name_tag] = L_B_1
        L_B_1_mean['%s' % get_name_tag] = L_B_1.mean()

        # Code for many different Bands:
        data_all_bands["%s" % get_name_tag] = []
        for k in Band_nrs[lowest_band:highest_band]: # Bands 8-11
            L_B = radiance_scales[k] * (Bands[k] - radiance_offsets[k])     # List with all bands
            print "Appending mean value of {} for band {} out of {}".format(L_B.mean(), Band_nrs[k], len(Band_nrs))
            data_all_bands['%s' % get_name_tag].append(L_B.mean())          # Mean radiance values

        i=i+1
        print "data added. Adding i+1 = ", i

    except AttributeError:
        print "\n*******************************"
        print "Can't open file {}".format(workdir1 + "/" + fe)
        print "Skipping this file..."
        print "*******************************"
        broken_files.append(workdir1 + "/" + fe)
        i=i+1

Answer 1

在不知道您的确切数据源和所需输出等的情况下，很难给您一个具体的答案。话虽如此，您似乎拥有 MODIS 图像的原生 .hdf 格式，并希望进行一些子集化以将图像引用到同一区域，然后进行绘图等。

查看gdal 模块中的gdal.Warp() 可能会对您有所帮助。此方法能够获取.hdf 文件并将一系列图像子集到具有相同分辨率/行数和列数的相同边界框。

然后您可以分析和绘制这些图像/比较像素等。

我希望这可以为您提供一个良好的入门起点。

gdal.Warp 文档：https://gdal.org/python/osgeo.gdal-module.html#Warp

更一般的warp帮助：https://www.gdal.org/gdalwarp.html

类似这样的：

import gdal

# Set up the gdal.Warp options such as desired spatial resolution,
# resampling algorithm to use and output format.
# See: https://gdal.org/python/osgeo.gdal-module.html#WarpOptions
# for other options that can be specified.
warp_options = gdal.WarpOptions(format="GTiff",
                                outputBounds=[min_x, min_y, max_x, max_y],
                                xRes=res,
                                yRes=res,
                                # PROBABLY NEED TO SET t_srs TOO
                                )

# Apply the warp.
# (output_file, input_file, options)
gdal.Warp("/path/to/output_file.tif",
          "/path/to/input_file.hdf",
          options=warp_options)

要编写的确切代码：

# Apply the warp.
# (output_file, input_file, options)
gdal.Warp('/path/to/output_file.tif',
          '/path/to/HDF4_EOS:EOS_SWATH:"MYD021KM.A2009034.1345.006.2012058160107.hdf":MODIS_SWATH_Type_L1B:EV_1KM_RefSB',
          options=warp_options)