Работа с растровыми данными

Работа с геопространственными данными в Python

Joris Van den Bossche

Open source software developer and teacher, GeoPandas maintainer

Иллюстрация сетки значений

Источник изображения: документация QGIS

Растровые данные

Карта вероятности осадков

Растровые данные с несколькими каналами

Пакет rasterio

import rasterio

«Pythonic»-обёртка над GDAL
Чтение и запись растровых файлов
Инструменты обработки (маскирование, перепроецирование, ресемплинг, ..)

https://rasterio.readthedocs.io/en/latest/

Открытие растрового файла

import rasterio

src = rasterio.open("DEM_world.tif")

Метаданные:

src.count

src.width, src.height

(4320, 2160)

Растровые данные = массив numpy

array = src.read()

Стандартный массив numpy:

array

array([[[-4290, -4290, -4290, ..., -4290, -4290, -4290],
        [-4278, -4278, -4278, ..., -4278, -4278, -4278],
        [-4269, -4269, -4269, ..., -4269, -4269, -4269],
        ...,
        [ 2804,  2804,  2804, ...,  2804,  2804,  2804],
        [ 2804,  2804,  2804, ...,  2804,  2804,  2804],
        [ 2804,  2804,  2804, ...,  2804,  2804,  2804]]], dtype=int16)

Визуализация растрового набора данных

Метод rasterio.plot.show():

import rasterio.plot

rasterio.plot.show(src, cmap='terrain')

Извлечение информации на основе векторных данных

rasterstats: сводная статистика растровых геопространственных данных на основе векторных геометрий (https://github.com/perrygeo/python-rasterstats)

Извлечение растровых значений с помощью rasterstats

Для точечных векторов:

rasterstats.point_query(geometries, "path/to/raster", 
                        interpolation='nearest'|'bilinear')

Для полигональных векторов:

rasterstats.zonal_stats(geometries, "path/to/raster",
                        stats=['min', 'mean', 'max'])

Извлечение растровых значений с помощью rasterstats

result = rasterstats.zonal_stats(countries.geometry, "DEM_gworld.tif", 
                                 stats=['mean'])

countries['mean_elevation'] = pd.DataFrame(result)

countries.sort_values('mean_elevation', ascending=False).head()

            name    continent                     geometry  mean_elevation
157   Tajikistan         Asia  POLYGON ((74.98 37.41, ...      3103.231105
85    Kyrgyzstan         Asia  POLYGON ((80.25 42.34, ...      2867.717142
24        Bhutan         Asia  POLYGON ((91.69 27.77, ...      2573.559846
119        Nepal         Asia  POLYGON ((81.11 30.18, ...      2408.907816
6     Antarctica   Antarctica  (POLYGON ((-59.57 -80.04...     2374.075028
..           ...          ...                          ...             ...

Давайте потренируемся!

Работа с геопространственными данными в Python