PostGIS для хранения данных на решетке

[kaa]

Спасибо, звучит похоже на то что надо. А какой самый разумный способ сгенерить растр (GeoTIFF) для таких данных ? Есть уже готовые утилиты или надо самому реализовывать разбор данных и их конвертацию?

[Ariki]

Нужно использовать библиотеку GDAL или утилиты, входящие в её поставку. Например, сформировать виртуальный набор данных, чтобы получить доступ к сырым бинарным данным, и использовать утилиту gdal_translate для конвертации его в GeoTIFF с нужными параметрами.

[SergeyRyzhkov]

Ariki
Кстати тоже вариант. Присоединяюсь к "спасибо" !

А Вы могли бы привести пример как сделать тифф по вышеприведенным данным?
И еще вопрос по данному решению, в чем преимущество перед просто хранить в базе?
По размеру не вижу выигрыш, по возможности рассчитывать в режиме online (если нужно будет считать что-то другое) - только минус.

[Ariki]

Тут преимущество в том, что не нужно читать весь растр целиком, то есть обработчику запроса HTTP не требуется выделять много памяти. Вдобавок решение получается дружелюбным к кэшированию: если скорость доступа критична, можно воткнуть какой-нибудь Redis для временного хранения часто используемых тайлов.
Если вы имеете в виду преимущества растровой формы представления над таблицей значений в БД, то она гораздо эффективней, так как не требует хранения координат каждого значения и к тому же хорошо сжимается алгоритмами компрессии.
Но эффективность, конечно, нужно проверять на практике.
Пример конвертации в TIFF постараюсь сделать вечером.

[kaa]

попробовал сделать tif вот такой командой:

Код: Выделить всё

gdal_translate cmorph.vrt cmorph.tif

с вот таким c vrt файлом

Код: Выделить всё

<VRTDataset rasterXSize="4948" rasterYSize="1649">
  <VRTRasterBand dataType="CFloat32" band="1" subClass="VRTRawRasterBand">
    <SourceFilename relativetoVRT="1">cmorph/CMORPH_V0.x_RAW_8km-30min_2015050412</SourceFilename>
    <ImageOffset>0</ImageOffset>
    <PixelOffset>4</PixelOffset>
    <LineOffset>19792</LineOffset>
    <ByteOrder>LSB</ByteOrder>
  </VRTRasterBand>
</VRTDataset>

вопрос а как тогда к этому файлу привязать геоинформацию? gdal_edit похоже требует tfw файл. Как его можно софрмировать для таких данных?

[Ariki]

Эта статья должна вам помочь.
Кстати, эти коэффициенты можно указать прямо в файле vrt, в элементе GeoTransform.
Ещё можно воспользоваться утилитой gdalwarp с параметром -te

[Ariki]

Чтобы оптимизировать GeoTIFF, добавьте опции к gdal_translate:

Код: Выделить всё

gdal_translate -co TILED=YES -co INTERLEAVE=BAND -co COMPRESS=LZW cmorph.vrt cmorph.tif

[kaa]

Долгота меняется от 0.036378335 с шагом 0.072756669, широта от -59.963614 с шагом 0.072771377.
Соответственно задаю преобразование в vrt

Код: Выделить всё

<GeoTransform>0.036378335,  0.072756669,  0.0,  -59.963614,  0.0, 0.072771377</GeoTransform>

Верно? Но как в таком случае быть с тем что что долгота после 180 должна стать отрицательной? И нужно ли в растр как то добавлять проекцию?

Получившийся растр не добавляется в PostGIS:

Код: Выделить всё

 raster2pgsql -f rast cmorph.tif precip_cmorph_raster
Processing 1/1: cmorph.tif
ERROR: convert_raster: The pixel type of band 1 is a complex data type.  PostGIS raster does not support complex data types
ERROR: process_rasters: Could not process raster: cmorph.tif
ERROR: Unable to process rasters

[Ariki]

Зависит от того, как вы потом будете обращаться к этому растру. Скорее всего, лучшим вариантом будет прописать ему соответствующую проекцию, а потом перепроецировать, чтобы нулевой меридиан оказался посередине. Как будет выглядеть проекция - сходу не напишу, надо поэкспериментировать.

[Ariki]

A gdalinfo cmorph.tif что показывает?

[kaa]

Код: Выделить всё

gdalinfo cmorph.tif
Driver: GTiff/GeoTIFF
Files: cmorph.tif
Size is 4948, 1649
Coordinate System is `'
Origin = (0.036378335000000,-59.963614000000000)
Pixel Size = (0.072756669000000,0.072771377000000)
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (   0.0363783, -59.9636140)
Lower Left  (   0.0363783,  60.0363867)
Upper Right (     360.036,     -59.964)
Lower Right (     360.036,      60.036)
Center      ( 180.0363774,   0.0363863)
Band 1 Block=256x256 Type=CFloat32, ColorInterp=Gray

[Ariki]

Что-то я туплю. У вас же там и правда комплексные значения пикселей (CFloat32). Если так и должно быть, то, видимо, PostGIS Raster вам не подходит. Или там на самом деле должно быть два канала с вещественными значениями? Есть где-нибудь документация по формату исходных файлов?

[kaa]

Извиняюсь, это я туплю. Конечно же там нужен был Float32. Поменял, стал создаваться нормальный tif.

Зависит от того, как вы потом будете обращаться к этому растру. Скорее всего, лучшим вариантом будет прописать ему соответствующую проекцию, а потом перепроецировать, чтобы нулевой меридиан оказался посередине. Как будет выглядеть проекция - сходу не напишу, надо поэкспериментировать.

В идеале хотелось бы чтобы можно было одним запросом считать среднее значение по территории, что то вроде:

Код: Выделить всё

select ST_SummaryStats(ST_Clip(raster, ST_GeomFromText('POLYGON(( ... ))')), true);

Где территория задается реальными координатами.

[Ariki]

Ух ты, я и не знал, что PostGIS Raster умеет считать статистику. Действительно, это лучше всего предложенного выше.
Может быть, эта ссылка вам поможет перенести начальный меридиан.

[kaa]

Спасибо огромное.
Переносить начальный меридиан не потребовалось, оказалось достаточно оборачивать геометрию в ST_Shift_Longitude (она конвертирует долготу из представления -180..180 в 0..360). В общем растр в PostGIS лег, занимает всего 1мб (может из за компрессии). Буду изучать насколько получающиеся значения адекватны.