MODIS AQUA Data - Укладка данных мозаики с помощью Python GDAL

Question

Я знаю, как получить доступ и построить поднаборы данных, используя gdal и python.Однако мне интересно, есть ли способ использовать данные GEO, содержащиеся в файле HDF4, чтобы я мог смотреть на одну и ту же область в течение многих лет.

И, если возможно, можно ли вырезать область данных и как это сделать?

ОБНОВЛЕНИЕ:

Если быть более точным: я построил данные MODIS и, как вы можете видеть,ниже река движется вниз (прямоугольное строение, верхний левый угол).Таким образом, в течение всего года это не то же самое место, которое я наблюдаю.

В подкатодах есть каталог, называемый Полями геолокации с Long и Altкаталоги.Так можно ли получить доступ к этой информации или наложить ее на данные, чтобы вырезать определенную область?

Если мы, например, взглянем на изображение НАСА ниже, можно ли будет сократить его между 10-15чередующийсяи длиной от -5 до 0.

Вы можете загрузить образец файла, скопировав URL ниже:

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

Не зная вашего точного источника данных, желаемого результата и т. Д., Трудно дать вам конкретный ответ.С учетом вышесказанного создается впечатление, что у вас есть собственный формат .hdf изображений MODIS и вы хотите выполнить некоторые поднаборы, чтобы изображения были привязаны к той же области, затем построить график и т. Д.

Это может помочь вампосмотрите на gdal.Warp() из модуля gdal.Этот метод может взять файл .hdf и поместить серию изображений в один и тот же ограничивающий прямоугольник с одинаковым разрешением / количеством строк и столбцов.

Затем вы можете проанализировать и построить эти изображения / сравнить пиксели и т. Д.

Я надеюсь, что это даст вам хорошую отправную точку для начала.

gdal.Warp документы: https://gdal.org/python/osgeo.gdal-module.html#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)