Как я могу открыть определенный набор данных (полосу) из изображения MODIS из нескольких наборов данных с помощью Rasterio? Rasterio IndexError: индекс полосы вне диапазона

Как открыть определенный набор данных из изображения MODIS из нескольких наборов данных с помощью Rasterio?

Я разместил некоторые примеры данных на GitHub: https://github.com/SteveObert/rasterIO_question/tree/master/data

Если я открою файл MODIS HDF только с одной полосой, приведенный ниже код будет работать так, как я хочу:

import rasterio

rasterfileMulti = 'MOD10A1_multiband_HEGOUT.hdf'
rasterfileSingle = 'MOD10A1_singleband_HEGOUT.hdf'
shapeFile = 'SFP_drainage.shp'


# Read an HDF into an array

dataset = rasterio.open('rasterfileSingle')
band1 = dataset.read(1)
print(band1)

Выход:

array([[ 53,  53, 250, ..., 255, 255, 255],
   [ 56,  56,  56, ..., 255, 255, 255],
   [ 56,  56,  49, ..., 255, 255, 255],
   ...,
   [ 78,  78,  78, ...,  53,  50,  50],
   [ 72,  78,  78, ...,  57,  57,  57],
   [ 72,  72,  72, ...,  61,  61,  61]], dtype=uint8)

Однако, если я пытаюсь открыть файл MODIS HDF с несколькими наборами данных, я получаю ошибку "Rasterio IndexError: диапазонный диапазон вне диапазона", как показано ниже.

rasterfileMulti = MOD10A1_multiband_HEGOUT.hdf

dataset2 = rasterio.open('rasterfileMulti')
band1 = dataset.read(1)
print(band1)

Ошибка из приведенного выше кода выглядит следующим образом:

/Users/steve/anaconda3/lib/python3.6/site-packages/rasterio/__init__.py:193: UserWarning: Dataset has no geotransform set.  Default transform will be applied (Affine.identity())
s.start()
Traceback (most recent call last):

File "<ipython-input-9-584312f89d76>", line 3, in <module>
band1 = dataset.read(1)

File "rasterio/_io.pyx", line 720, in rasterio._io.RasterReader.read

IndexError: band index out of range

В конечном итоге я хочу обрезать растр в шейп-файл. Приведенный ниже код работает так, как я хочу, если изображение Modis имеет только одну полосу, в данном случае NDSI_Snow_Cover.

import fiona
import rasterio
rasterfileSingle = MOD10A1_singleband_HEGOUT.hdf
shapeFile = SFP_drainage.shp

with fiona.open(shapeFile, 'r') as shapefile:
    features = [feature['geometry'] for feature in shapefile]


with rasterio.open(rasterfileMulti) as src:
    out_image, out_transform = rasterio.mask.mask(src, features,
                                                        crop=True)
    out_meta = src.meta.copy()


out_meta.update({'driver': 'GTiff',
                 'height': out_image.shape[1],
                 'width': out_image.shape[2],
                 'transform': out_transform})
with      rasterio.open('/Users/steve/Documents/classes/Geos_505/project_Payette/working/data_files/test_clip_out.tif', 'w', **out_meta) as dest:
    dest.write(out_image)
Источник

1 ответ

Решение

Во-первых, нам нужно знать имя набора данных. gdalinfo предоставляет следующие 7 поднаборов:

Subdatasets:
  SUBDATASET_1_NAME=HDF4_EOS:EOS_GRID:"MOD10A1_multiband_HEGOUT.hdf":MOD_Grid_Snow_500m_16:NDSI_Snow_Cover
  SUBDATASET_1_DESC=[350x831] NDSI_Snow_Cover MOD_Grid_Snow_500m_16 (8-bit unsigned integer)           
  SUBDATASET_2_NAME=HDF4_EOS:EOS_GRID:"MOD10A1_multiband_HEGOUT.hdf":MOD_Grid_Snow_500m_16:NDSI_Snow_Cover_Basic_QA
  SUBDATASET_2_DESC=[350x831] NDSI_Snow_Cover_Basic_QA MOD_Grid_Snow_500m_16 (8-bit unsigned integer)
  SUBDATASET_3_NAME=HDF4_EOS:EOS_GRID:"MOD10A1_multiband_HEGOUT.hdf":MOD_Grid_Snow_500m_16:NDSI_Snow_Cover_Algorithm_Flags_QA
  SUBDATASET_3_DESC=[350x831] NDSI_Snow_Cover_Algorithm_Flags_QA MOD_Grid_Snow_500m_16 (8-bit unsigned integer)
  SUBDATASET_4_NAME=HDF4_EOS:EOS_GRID:"MOD10A1_multiband_HEGOUT.hdf":MOD_Grid_Snow_500m_16:NDSI
  SUBDATASET_4_DESC=[350x831] NDSI MOD_Grid_Snow_500m_16 (16-bit integer)
  SUBDATASET_5_NAME=HDF4_EOS:EOS_GRID:"MOD10A1_multiband_HEGOUT.hdf":MOD_Grid_Snow_500m_16:Snow_Albedo_Daily_Tile
  SUBDATASET_5_DESC=[350x831] Snow_Albedo_Daily_Tile MOD_Grid_Snow_500m_16 (8-bit unsigned integer)
  SUBDATASET_6_NAME=HDF4_EOS:EOS_GRID:"MOD10A1_multiband_HEGOUT.hdf":MOD_Grid_Snow_500m_16:orbit_pnt
  SUBDATASET_6_DESC=[350x831] orbit_pnt MOD_Grid_Snow_500m_16 (8-bit integer)
  SUBDATASET_7_NAME=HDF4_EOS:EOS_GRID:"MOD10A1_multiband_HEGOUT.hdf":MOD_Grid_Snow_500m_16:granule_pnt
  SUBDATASET_7_DESC=[350x831] granule_pnt MOD_Grid_Snow_500m_16 (8-bit unsigned integer)

Теперь вы можете открыть любой поднабор, передав в качестве аргумента полное имя GDAL поднабора. Например:

rasterio.open('HDF4_EOS:EOS_GRID:"MOD10A1_multiband_HEGOUT.hdf":MOD_Grid_Snow_500m_16:NDSI_Snow_Cover')

Изменить: я добавил внешние одинарные кавычки в соответствии с предложением авторов.

Источник
Другие вопросы по тегам