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

Эта страница опубликована в основном списке статей сайта
по адресу http://gis-lab.info/qa/postgis-rasters-example.html

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

ftp://ftp.cpc.ncep.noaa.gov/precip/CMORPH_V0.x/RAW/8km-30min/

. Это данные об осадках рассчитанных по алгоритму CMORPH. Файл с данными на сервере содержит 2 временных среза для каждого часа начиная с октября 2014 (за первые полчаса и вторые полчаса). Данные одного среза это массив 4948 x 1649 вещественных чисел рассчитанных для сетки с шагом 0.072756669 по долготе и 0.072771377 по широте начиная от 0.036378335E, -59.963614S.

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

Понятно, что все эти задачи можно решать напрямую, читая данные из файлов, суммируя значения. Но возник вопрос: а может ли для этой задачи помочь PostGIS? Можно ли как то занести подобные данные в базу, чтобы потом решать поставленные задачи простыми вызовами PostGIS SQL?

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

CREATE TABLE cmorph_30min_8km
(
  id bigint NOT NULL,
  measured_on timestamp with time zone NOT NULL,
  value real NOT NULL,
  location geometry(Point,4326) NOT NULL,
  CONSTRAINT cmorph_30min_8km_pkey PRIMARY KEY (id)
);
CREATE INDEX gist_cmorph_30min_8km
  ON cmorph_30min_8km
  USING gist
  (location);

Один временной cрез отображается в ~10M строк в PostGIS, что приводит к чудовищному размеру базы (около 1gb дискового пространства на срез, из которых половину занимает индекс). Попытка вместо геометрии хранить просто положение точки в решетке(int x и int y) ситуацию не улучшила, размер данных получался такой же.

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

1. Отдельная таблица определяющую сетку (узлы).
   • Идентификатор
   • Геометрия точки
   • Пространственный индекс естественно
2. Отдельная таблица для значений
   • Идентификатор
   • Внешний ключ на табл 1
   • Дата измерения
   • Значение

Это ситуацию улучшило, но не кардинально. Таблица с координатами занимала все равно 1gb, но каждый срез обходился уже всего в 400mb.

Возникла идея использовать растровый формат для хранения этих данный. Растр из исходного файла можно получить с помощью утилиты gdal_translate. Эта утилита может использовать бинарный файл как источник данных для растра. То, каким способом интерпретируются данные - можно описать с помощью GDAL Virtual Format. Для этого формируются два VRT файла - cmorph00.vrt и cmorph30.vrt (так как каждый исходный файл содержит два среза, за первые и вторые полчаса) по шаблону:

<!-- задается размер получаемого растра 4948 на 1649 -->
<VRTDataset rasterXSize="4948" rasterYSize="1649">
  <!-- коэффициенты афинного преобразования из положения пиксела (колонка, строка) в геокоординаты. Подробнее здесь: http://www.gdal.org/gdal_datamodel.html -->
  <GeoTransform>0.036378335,  0.072756669,  0.0,  -59.963614,  0.0, 0.072771377</GeoTransform>
  <!-- тип данных и номер слоя в получаемом растре -->
  <VRTRasterBand dataType="Float32" band="1" subClass="VRTRawRasterBand">
    <!-- значение интерпретируемое как отсутствие данных -->
    <NoDataValue>-999.0</NoDataValue>
    <!-- исходный файл -->
    <SourceFilename relativetoVRT="1">cmorph/CMORPH_V0.x_RAW_8km-30min_2015050519</SourceFilename>
    <!-- смещение в байтах относительно начала файла (0 для первого среза, 4948*1649*4 = 32637008 для второго) -->
    <ImageOffset>0</ImageOffset>
    <!-- смещение для каждого пиксела  -->
    <PixelOffset>4</PixelOffset>
    <!-- смещение для каждой линии: 4948*4 = 19792 -->
    <LineOffset>19792</LineOffset>
    <!-- формат бинарного представления чисел: LSB - младшие байты первые, MSB - старшие байты первые -->
    <ByteOrder>LSB</ByteOrder>
  </VRTRasterBand>
</VRTDataset>

Затем из исходного файла формируются два файла GeoTIFF размером всего по 600kb:

gdal_translate -co INTERLEAVE=BAND -co COMPRESS=LZW cmorph00.vrt cmorph00.tif

gdal_translate -co INTERLEAVE=BAND -co COMPRESS=LZW cmorph30.vrt cmorph30.tif

Параметры INERLEAVE и COMPRESS требуются для создания сжатого и оптимального для загрузки tif файла.

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

gdal_calc.py --co="COMPRESS=LZW" -A cmorph1.tif -B cmorph2.tif --outfile=result.tif --calc="(A+B)/2"

Деление на 2 тут требуется так как исходные данные имеют размерность мм/час, но взяты за 30 мин.

Полученный растр можно загрузить в PostGIS, подготовив SQL вызовом

raster2pgsql -N -999.0 -f rast cmorph00.tif cmorph_8km_30min > cmorph.sql

. В том числе есть возможность создать скрипт без предварительного создания растра, прямо из VRT файла:

raster2pgsql -N -999.0 -f rast cmorph00.vrt cmorph_8km_30min > cmorph.sql

После выполнения полученного cmorph.sql в базе создается таблица cmorph_8km_30min содержащая нужный растр. В PostGIS, то он будет занимать менее 1mb дискового пространства, что в 1000 раз лучше предыдущих подходов. Но что можно делать с данными в растровом представлении кроме как показывать их как картинку? Оказывается PostGIS включает встроенные аналитические функции для работы с растрами. Например функция ST_Histogram - возвращает распределение значений по интервалам, ST_SummaryStats позволит выбрать минимальное, максимальное и среднее значение для заданного растра. То есть для решения изначальной задачи о количестве осадков в определенной территории в определенном временном интервале достаточно выполнить следующий запрос к базе:

select (stats).mean as mean, measured_on from (select ST_SummaryStats(ST_Clip(rast, ST_Shift_Longitude(ST_GeomFromText('POLYGON(( ... ))'))), true) as stats, measured_on from cmorph_8km_30min where measured_on >= ... and measured_on <= ...) q

Здесь POLYGON(( ... )) задает требуемую территорию, а ST_Shift_Longitude используется для конверсии долготы с интервала -180..180 на 0..360.

Полный список функций PostGIS для работы с растрами можно посмотреть здесь: http://postgis.net/docs/RT_reference.html .