Растровый формат данных - Геоаналитика на Python

Растровые пространственные данные — это способ представления географической информации в виде регулярной сетки (матрицы) одинаковых по размеру ячеек — пикселей. Каждая ячейка соответствует определённой территории и содержит числовое значение (или набор значений), например: высоту, температуру, интенсивность отражённого излучения и др.

Растровый формат особенно удобен для данных, которые непрерывно изменяются в пространстве, то есть описывают различные поверхности: рельеф, осадки, температуру, влажность почвы и др.

Размер пикселя определяет пространственное разрешение данных: чем он меньше, тем более детально растр описывает территорию.

С растровым форматом пространственных данных мы уже кратко познакомились в первом разделе.

В этом разделе мы поработаем с реальными данными о плотности населения с портала WorldPop

0. Импорт библиотек¶

import rasterio
import numpy as np
import matplotlib.pyplot as plt

import osmnx as ox

1. Данные WorldPop¶

Данные о плотности населения на портале WorldPop представляют собой оценку распределения населения в виде растровой поверхности. Такие оценки получают, комбинируя данные переписей населения с различными вспомогательными источниками — например, спутниковыми снимками, данными о застройке и инфраструктуре.

В результате население «распределяется» по ячейкам растра в зависимости от характеристик территории: где выше вероятность проживания людей, там значения оказываются больше.

1.1 Экспортируем данные¶

Экспортируем данные WorldPop о распределении населения за 2020 год. Ниже приведён пример кода для скачивания растрового файла с официального сервера. Обратите внимание, что файл достаточно большой (около 300 МБ), поэтому загрузка может занять некоторое время.

Если вы не хотите тратить время на скачивание, готовый файл уже размещён в директории с материалами курса, и этот шаг можно пропустить.

В дальнейших примерах мы будем использовать уже подготовленный файл — обрезанный по территории Тульской области.

#Можно выгрузить данные с портала World Pop (300 Mb + )

# # Ссылка на данные WorldPop
# url = "https://data.worldpop.org/GIS/Population/Global_2000_2020_Constrained/2020/BSGM/RUS/rus_ppp_2020_constrained.tif"

# # Имя файла для сохранения
# output_file = "worldpop_russia_2020.tif"

# # Скачивание файла
# response = requests.get(url, stream=True)

# if response.status_code == 200:
#     with open(output_file, "wb") as file:
#         for chunk in response.iter_content(chunk_size=1024):
#             file.write(chunk)
#     print(f"File saved as {output_file}")
# else:
#     print("Failed to download file. Status code:", response.status_code)

1.2 Чтение растровых данных¶

Откроем подготовленный растр, обрезанный по границам Тульской области, с помощью библиотеки rasterio. Растр tula_region_population.tif можно найти в репозитории.

Значения в растре соответствуют оценке численности населения в пределах пикселя.

После открытия файла мы выведем его основные характеристики:

CRS (система координат) — показывает, в какой системе заданы данные;
границы (bounds) — пространственный охват растра;
разрешение (resolution) — размер одной ячейки (пикселя) в единицах системы координат.

output_file = './data/tula_region_population.tif'


with rasterio.open(output_file) as dataset:
    print("CRS:", dataset.crs)
    print("Bounds:", dataset.bounds)
    print("Resolution:", dataset.res)

CRS: EPSG:4326
Bounds: BoundingBox(left=35.89855540000002, bottom=52.94600317428952, right=38.95188872112002, top=54.85183649999999)
Resolution: (0.00083333333, 0.0008333333300002048)

1.3 Визуализация растра на карте¶

После чтения данных мы можем визуализировать растр. По сути, растровый слой — это обычная матрица значений, поэтому его можно отобразить как изображение с помощью библиотеки matplotlib.

В примере ниже мы отображаем первый (и в данном случае единственный) канал растра и добавляем цветовую шкалу, чтобы интерпретировать значения.

with rasterio.open(output_file) as dataset:
    data = dataset.read(1)  

    plt.figure(figsize=(10, 10))
    plt.imshow(data, cmap='viridis') 
    plt.colorbar(label="Population Density") 
    plt.title("Population Density Tula Region (WorldPop, Russia, 2020)")
    plt.show()

src.read(1) считывает первый канал растра в виде двумерного массива (матрицы);
plt.imshow(data, cmap='viridis') отображает матрицу как изображение, где значения пикселей преобразуются в цвета;
cmap='viridis' задаёт цветовую схему;
plt.colorbar() добавляет шкалу, которая показывает, каким значениям соответствуют цвета;
plt.title() добавляет заголовок для удобства интерпретации.

1.4 Обработка отсутствующих значений¶

На предыдущем шаге при визуализации растра можно было заметить, что изображение выглядит «искажённым»: значительная часть территории отображается в узком диапазоне цветов, а цветовая шкала смещена.

Это происходит потому, что в данных присутствует большое количество сильно отрицательных значений, которые влияют на автоматическое масштабирование цветовой шкалы.

Давайте разберёмся, что это за значения. Для этого посмотрим на базовую статистику растра:

with rasterio.open(output_file) as dataset:
    data = dataset.read(1) 
    # Выводим статистику
    print("Min value:", np.min(data))
    print("Max value:", np.max(data))
    print("Mean value:", np.mean(data))

Min value: -99999.0
Max value: 161.41058
Mean value: -97470.01

Мы увидим, что минимальное значение растра равно -99999. Очевидно, что это не может быть значением плотности населения.

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

Чтобы корректно работать с данными, такие пиксели нужно исключать из анализа.

1.4.1. Исключение NoData значений¶

Библиотека rasterio позволяет сразу учитывать отсутствующие значения с помощью masked array:

with rasterio.open(output_file) as dataset:
    data = dataset.read(1, masked=True)

В этом случае все пиксели с NoData автоматически маскируются и не участвуют в вычислениях и визуализации.

Прверим какое теперь минимальное значение у исследуемого растра

print("Min (masked):", data.min())

Min (masked): 0.66187334

1.4.2. Визуализация без NoData¶

Теперь посмотрим на растр, исключив отсутствующие значения:

plt.figure(figsize=(10, 10))
plt.imshow(data, cmap="viridis")
plt.colorbar(label="Population Density")
plt.title("Population Density Tula Region (WorldPop, Russia, 2020)")
plt.show()

2. Обрезка растра по векторному слою¶

В предыдущих шагах мы работали с растром, охватывающим всю Тульскую область. Однако часто требуется анализировать данные только для конкретной территории — например, отдельного города.

В таких случаях растровые данные обрезают (клиппируют) по границам интересующей области.

В качестве примера ограничим данные территорией города Тулы.

Для этого нам понадобятся границы города в виде векторного слоя. Мы загрузим их из OpenStreetMap с помощью библиотеки osmnx.

2.1. Экспорт векторных данных¶

Загрузим границу Тулы из OpenStreetMap

area_name = "Тула, Россия"

gdf_borders = ox.geocode_to_gdf(area_name)

gdf_borders.explore(tiles='cartodbpositron')

2.2 Обрезка растра по границам¶

Теперь обрежем растр по границе города Тулы.

Для этого будем использовать функцию mask из библиотеки rasterio, которая позволяет обрезать растровые данные по произвольной геометрии.

Перед обрезкой важно убедиться, что растр и векторные границы находятся в одной системе координат (CRS) — иначе результат будет некорректным.

2.2.1 Проверка систем координат¶

with rasterio.open(output_file) as dataset:
    
    # Если CRS различаются — перепроецируем границы
    if gdf_borders.crs != dataset.crs:
        gdf_borders = gdf_borders.to_crs(dataset.crs)

2.2.2. Подготовка геометрии¶

Функция mask принимает геометрию в формате GeoJSON, поэтому преобразуем границы:

import json

geometries = [json.loads(gdf_borders.to_json())['features'][0]['geometry']]

2.2.3. Обрезка растра¶

from rasterio.mask import mask

with rasterio.open(output_file) as dataset:
    out_image, out_transform = mask(dataset, geometries, crop=True)

out_image — обрезанный растр
out_transform — новое пространственное преобразование

Параметр crop=True означает, что растр будет уменьшен до границ геометрии.

2.2.4. Обновление метаданных¶

После обрезки необходимо обновить метаданные растра:

out_meta = dataset.meta.copy()
out_meta.update({
    "driver": "GTiff",
    "height": out_image.shape[1],
    "width": out_image.shape[2],
    "transform": out_transform
})

2.2.5. Сохранение результата¶

file_cropped_path = "./data/tula_cropped_population.tif"

with rasterio.open(file_cropped_path, "w", **out_meta) as dest:
    dest.write(out_image)

2.3. Анализ полученного растра¶

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

2.3.1 Основная информация о растре¶

with rasterio.open(file_cropped_path) as dataset:
    print("CRS:", dataset.crs)  
    print("Bounds:", dataset.bounds)  
    print("Width, Height:", dataset.width, dataset.height)  
    print("Number of bands:", dataset.count)  
    print("Data type:", dataset.dtypes)  
    print("Transform:", dataset.transform)

CRS: EPSG:4326
Bounds: BoundingBox(left=37.49522206028002, bottom=54.061836503159796, right=37.741888725960024, top=54.300169835539855)
Width, Height: 296 286
Number of bands: 1
Data type: ('float32',)
Transform: | 0.00, 0.00, 37.50|
| 0.00,-0.00, 54.30|
| 0.00, 0.00, 1.00|

Здесь мы проверяем:

CRS — систему координат растра;
Bounds — пространственные границы;
Width, Height — размер растра в пикселях;
Number of bands — количество каналов;
Data type — тип значений;
Transform — параметры геопривязки растра.

Обратите внимание, что растр находится в географической системе координат — EPSG:4326.

Это означает, что:

координаты заданы в градусах;
размер пикселя также выражен в градусах, а не в метрах.

Это важно учитывать при дальнейшем анализе, например при расчёте расстояний, площадей или плотности.

2.3.2 Размер пикселя¶

Размер одного пикселя можно получить напрямую с помощью атрибута res библиотеки rasterio:

with rasterio.open(file_cropped_path) as dataset:
    pixel_width, pixel_height = dataset.res
    
    print("Pixel Width:", pixel_width)
    print("Pixel Height:", pixel_height)

Pixel Width: 0.00083333333
Pixel Height: 0.0008333333300002048

pixel_width — размер пикселя по оси X;
pixel_height — размер пикселя по оси Y.

Важно помнить, что в данном случае размер пикселя задан в градусах, так как растр находится в системе координат EPSG:4326.

2.3.3 Визуализация обрезанного растра¶

Теперь посмотрим, как выглядит новый растровый слой после обрезки:

with rasterio.open(file_cropped_path) as dataset:
    data = dataset.read(1, masked=True)

plt.figure(figsize=(10, 10))
plt.imshow(data, cmap="viridis")
plt.colorbar(label="Population")
plt.title("Population Distribution — Tula City (WorldPop, 2020)")
plt.show()

3. Перепроецирование растра¶

Сейчас наш растровый слой находится в географической системе координат (EPSG:4326), где координаты и размер пикселя заданы в градусах.

Для большинства аналитических задач (расчёт расстояний, площадей, плотности) это неудобно, поэтому растр обычно перепроецируют в систему координат с метрическими единицами — например, в UTM.

3.1. Определение подходящей UTM-зоны¶

UTM-зона зависит от географического положения территории. Её можно определить автоматически на основе векторного слоя с границами.

Для этого воспользуемся возможностями библиотеки geopandas:

utm_crs = gdf_borders.estimate_utm_crs()
print("Estimated UTM CRS:", utm_crs)

Estimated UTM CRS: EPSG:32637

3.2. Перепроецирование растра¶

Для перепроецирования используем функции из библиотеки rasterio.

Перепроецирование выполняется отдельно для каждого канала растра, если каналов несколько

3.2.1. Подготовка параметров новой проекции¶

Сначала рассчитываем, как будет выглядеть растр в новой системе координат:

from rasterio.warp import calculate_default_transform

with rasterio.open(file_cropped_path) as dataset:
    transform, width, height = calculate_default_transform(
        dataset.crs, utm_crs, dataset.width, dataset.height, *dataset.bounds
    )

transform — новое аффинное преобразование
width, height — новые размеры растра

3.2.2. Обновление метаданных¶

Теперь обновим метаданные растра с учётом новой системы координат:

kwargs = dataset.meta.copy()
kwargs.update({
        "crs": utm_crs,
        "transform": transform,
        "width": width,
        "height": height
    })

3.2.3. Создание нового файла¶

Зададим путь для сохранения перепроецированного растра:

output_reprojected = "./data/tula_cropped_population_utm.tif"

3.2.4. Перепроецирование данных¶

Теперь выполняем перепроецирование и записываем результат в новый файл:

from rasterio.warp import reproject, Resampling

with rasterio.open(file_cropped_path) as src:
    with rasterio.open(output_reprojected, "w", **kwargs) as dst:
        for i in range(1, src.count + 1):
            reproject(
                source=rasterio.band(src, i),
                destination=rasterio.band(dst, i),
                src_transform=src.transform,
                src_crs=src.crs,
                dst_transform=transform,
                dst_crs=utm_crs,
                resampling=Resampling.nearest
            )

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

3.3. Проверка результата¶

with rasterio.open(output_reprojected) as dataset:
    print("CRS:", dataset.crs)
    print("Resolution:", dataset.res)

CRS: EPSG:32637
Resolution: (75.37731228269381, 75.37731228269381)

После перепроецирования:

система координат стала проекционной (UTM);
размер пикселя теперь задан в метрах;
растр готов для дальнейшего пространственного анализа.

Итог¶

В этом разделе мы познакомились с растровым форматом пространственных данных и рассмотрели основные этапы работы с ним на примере данных WorldPop.

Мы:

научились открывать и анализировать растровые файлы с помощью библиотеки rasterio;
выявили и обработали отсутствующие значения (NoData);
обрезали растр по границам интересующей территории;
перепроецировали данные.

В результате мы получили подготовленный растровый слой, который можно использовать для дальнейшего пространственного анализа — например, для оценки распределения населения, агрегации данных или расчёта показателей доступности.