Объединение нескольких растров GeoTIFF в один

Эта проблема возникает часто, когда вы хотите работать с большей областью, чем охватывает один файл, или вы извлекли некоторые спутниковые данные, но получили половину захваченного изображения для какой-то даты, а другую половину для другой даты. это означает, что для анализа всего региона вам пришлось бы делать это несколько раз (в нескольких частях региона), но зачем это делать, когда вы можете объединить эти растры или набор данных спутниковых изображений в один файл. Посмотрим как?

В этом уроке я покажу, как использовать утилиту Gdal для объединения нескольких плиток вместе. Я также объясню, как установить привязки Gdal и Gdal-python в вашей системе, чтобы беспрепятственно использовать утилиту.

Цель

Чтобы объединить несколько растровых файлов GeoTiff в один файл.

Требования

1. Привязки Gdal и Gdal-python
2. Ниже пакеты Python

• Glob — Чтобы перечислить все растровые файлы в каталоге
• Matplotlib — для построения растровых файлов/спутниковых изображений.
• Numpy — для чтения растровых данных в массив numpy.

Входные данные

Входными данными для этой задачи будет n количество растровых файлов, которые вы хотите объединить в один.

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

Я загрузил приведенные ниже растровые файлы для использования в качестве входных данных для этого урока и сохранил их в папку data/.

1. растровый файл 1 — s3://sentinel-cogs/sentinel-s2-l2a-cogs/43/R/EL/2023/1/S2A_43REL_20230102_0_L2A/B04.tif
2. растровый файл 2 — s3://sentinel-cogs/sentinel-s2-l2a-cogs/43/R/EL/2023/1/S2B_43REL_20230104_0_L2A/B04.tif

Два файла были захвачены в две близкие даты 2023–01–02 и 2023–01–04, предоставляя частичные данные дозорной плитки для каждой даты. Эти частичные наборы данных объединяются, чтобы получить набор данных для полной дозорной плитки.

Установка и настройка Gdal

Gdal — это библиотека-транслятор для растровых и векторных форматов геопространственных данных.

Установка Gdal отличается и сложнее, чем установка обычного пакета/модуля Python. Это также зависит от операционной системы, которую вы настроили на своем компьютере, давайте установим gdal и его привязки на разные ОС:

Установка Gdal в Ubuntu:

Откройте терминал на своем компьютере и выполните следующие действия:

Чтобы получить последнюю версию GDAL, добавьте PPA в свои исходники, затем установите пакет gdal-bin (это должно автоматически получить все необходимые зависимости, включая, по крайней мере, соответствующую версию libgdal).

sudo add-apt-repository ppa:ubuntugis/ppa

Как только вы добавите репозиторий, обновите исходные пакеты.

sudo apt-get update

Теперь вы сможете установить пакет GDAL.

sudo apt-get install gdal-bin

Чтобы проверить установку, вы можете запустить ogrinfo --version. Вам потребуется версия GDAL для установки правильных привязок Python.

ogrinfo --version
GDAL 3.6.2, released 2023/01/02

Перед установкой библиотек GDAL Python вам необходимо установить библиотеки разработки GDAL.

sudo apt-get install libgdal-dev

Вам также потребуется экспортировать пару переменных окружения для компилятора.

export CPLUS_INCLUDE_PATH=/usr/include/gdal
export C_INCLUDE_PATH=/usr/include/gdal

Теперь вы можете использовать pip для установки привязок Python GDAL.

pip install GDAL==<GDAL VERSION FROM ogrinfo --version>

Установка Gdal на Mac-OS:

Установите Homebrew, если у вас его нет. Скопируйте и вставьте следующие строки в свой терминал, чтобы установить доморощенный.

sr/bin/ruby -e "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/master/install)"

Сначала установите файлы заголовков, чтобы избежать сбоя GDAL.

brew install gdal --HEAD

Теперь установите GDAL

brew install gdal

Теперь вы можете использовать pip для установки привязок Python GDAL.

pip install GDAL==<GDAL VERSION FROM ogrinfo --version>

Проверить набор данных

Я покажу вам, как выглядят эти два входных растровых файла, визуализируя оба набора данных после их чтения в 2D-массив. Я собираюсь использовать приведенный ниже скрипт для чтения растрового файла GeoTiff в массив.

def rater_to_array(raster_file, output_data_type):
    ds = gdal.Open(geotif_file)
    number_of_bands = ds.RasterCount
    if number_of_bands == 1:
        raster_band = ds.GetRasterBand(1)
        no_data_value = raster_band.GetNoDataValue()
        array = band.ReadAsArray()
        array[np.where(array==)]=np.nan
        array = array.astype(output_data_type)
    else:
        array_rows = ds.RasterYSize
        array_cols = ds.RasterXSize
        array = np.zeros((array_rows,array_cols,number_of_bands))
        for i in range(1, number_of_bands+1):
            raster_array = None
            raster_band = ds.GetRasterBand(1)
            raster_array = band.ReadAsArray()
            raster_array[np.where(raster_array==no_data_value)]=np.nan
            array[...,i-1] = raster_array
    return array

Давайте визуализируем оба растровых файла и проверим, как это выглядит.

from osgeo import gdal
import numpy as np
import matplotlib.pyplot as plt
raster1 = 'data/S2A_43REL_20230102_0_L2A/B04.tif'
raster2 = 'data/S2B_43REL_20230104_0_L2A/B04.tif'
raster_arr1 = raster_to_array(raster1, np.float)
plt.subplot(2,2,1)
plt.imshow(raster_arr1)
raster_arr2 = raster_to_array(raster2, np.float)
plt.imshow(raster_arr2)
plt.set_cmap(colormap); 
plt.show()

Как вы видите два вышеприведенных входных растровых файла, эти спутниковые снимки сделаны в две разные даты, которые являются не чем иным, как двумя частями одного тайла. Я покажу вам, как мы можем использовать Gdal для объединения этих двух частей в один набор данных.

Сначала сохраните все файлы, которые вы хотите объединить, в список, используя glob. распечатайте список и посмотрите.

import glob
all_raster_files = glob.glob('data/*/*.tif')
print(all_raster_files)
['data/S2A_43REL_20230102_0_L2A/B04.tif',
 'data/S2B_43REL_20230104_0_L2A/B04.tif']

Теперь, когда у меня есть список файлов, я собираюсь написать ниже функцию для объединения этих файлов вместе и преобразования в один растр, используя утилиты gdal Warp, BuildVrt и Translate.

def merge_rasters(raster_file_list, width=None, height=None): ## Adding width and height as custom parameters if want to change the size of raster
      output_file = os.path.join('data', 'merged.tif')
      ds_lst = list()
      for raster in raster_file_list:
          ds = gdal.Warp('', raster, format='vrt', dstNodata=0,
                       dstSRS="+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0",
                       width=out_width, height=out_height)
          ds_lst.append(ds)
          del ds
      dataset = gdal.BuildVRT('', ds_lst, VRTNodata=0, srcNodata=0)
      ds1 = gdal.Translate(output_file, dataset)
      del ds1  
      del dataset
      return output_file

Сначала мы сохранили набор данных из каждого растрового файла в формате vrt, а затем использовали BuildVRT для объединения этого списка наборов растровых данных в один файл vrt. Далее использовалась утилита перевода для преобразования этого окончательного набора данных в формат GeoTiff.

Теперь я собираюсь визуализировать окончательный выходной файл data/merged.tif с помощью matplotlib.

from osgeo import gdal
import numpy as np
import matplotlib.pyplot as plt
merged_raster_path = 'data/merged.tif'
raster_arr = raster_to_array(merged_raster_path, np.float)
plt.imshow(raster_arr)
plt.show()

Давайте посмотрим их все вместе!

Большой ! Похоже, мы добились своих результатов. Как видите, я объединил файлы в файл merged.tif. Здесь следует отметить, что Гдаль работал здесь очень эффективно и мощно. Я бы посоветовал вам прочитать другие статьи, в которых у меня есть опыт использования утилит Gdal для решения различных проблем, связанных с набором данных растровых/GeoTiff или спутниковых изображений.