Пересечение растра с полигоном с помощью PostGIS - ошибка артефакта

15

Я использую PostGIS2.0 для пересечения растров и полигонов. Мне трудно понять, какую операцию мне следует использовать, и какой это самый быстрый способ выполнить. Моя проблема заключается в следующем:

  • У меня есть многоугольник и растр
  • Я хочу найти все пиксели, которые попадают в многоугольник, и получить сумму значений пикселей
  • И (обновленная проблема): я получаю массивные значения для некоторых пикселей, которых нет в исходном растре, когда я выполняю запрос

Мне трудно понять, следует ли мне использовать ST_Intersects()или ST_Intersection(). Я также не знаю, как лучше всего суммировать мои пиксели. Вот первый подход, который я попробовал (# 1):

SELECT 
    r.rast 
FROM
    raster as r, 
    polygon as p
WHERE
    ST_Intersects(r.rast, p.geom)

Это возвращает список rastзначений, с которыми я не уверен, что делать. Я попытался вычислить итоговую статистику, используя, ST_SummaryStats()но я не уверен, является ли это взвешенной суммой всех пикселей, которые лежат в пределах многоугольника.

SELECT  
        (result).count,
        (result).sum    
FROM (  
         SELECT 
            ST_SummaryStats(r.rast) As result
         FROM
            raster As r, 
            polygon As p
         WHERE
            ST_Intersects(r.rast, p.geom)    
    ) As tmp

Другой подход, который я попробовал (# 2), использует ST_Intersection():

 SELECT
        (gv).geom,         
        (gv).val
 FROM 
 (
    SELECT 
        ST_Intersection(r.rast, p.geom) AS gv
    FROM 
        raster as r, 
        polygon as p           
    WHERE 
        ST_Intersects(r.rast, p.geom)

      ) as foo;

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

Мне непонятно, какой порядок действий также самый быстрый. Должен ли я всегда выбирать raster, polygonили polygon, raster, или конвертировать полигон в растр, чтобы он был raster, raster?

РЕДАКТИРОВАТЬ: Я обновил подход № 2 с некоторыми подробностями из R.K.ссылки.

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

введите описание изображения здесь

И вот результат пересечения с использованием PostGIS:

введите описание изображения здесь

Проблема с результатом заключается в том, что возвращаются значения 21474836, которых нет в исходном растре. Я понятия не имею, почему это происходит. Я подозреваю, что это где-то связано с небольшими числами (деление почти на 0), но возвращает неправильный результат.

djq
источник
Что касается второго пункта, вы хотите получить сумму значений пикселей, которые пересекают многоугольник?
RK
Да, я использовал ST_SummaryStats()для № 1, но не уверен, как это сделать для № 2.
djq
Разместил ссылку на ссылку. Надеюсь, это поможет.
РК
2
Максимальное значение масштаба на вашей карте - это максимум 32-разрядного целого числа со знаком. Я не знаю, что это значит в вашем случае, но это может быть связано со значениями NA. Диапазон в вашем запросе может иметь нулевые значения, которые не обрабатываются должным образом. en.wikipedia.org/wiki/2147483647#2147483647_in_computing
yellowcap
6
FWIW, 21474836 - не максимальное значение 32-битного целого со знаком . Тем не менее, 2 ^ 31-1 = 2147483647 является максимумом, и обратите внимание, что 21474836 = 2147483647/100 (целочисленное деление). Это намекает на то, что внутренне некоторые NA создаются (возможно, вдоль граничных ячеек), они представляются как 2 ^ 31-1, а затем код «забывает», что это NA и (возможно, в процессе повторной выборки?) Делит их на 100.
whuber

Ответы:

6

Я нашел учебник по пересекающимся векторным полигонам с большим растровым покрытием, используя PostGIS WKT Raster . Может иметь ответы, которые вы ищете.

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

 CREATE TABLE caribou_srtm_inter AS
 SELECT id, 
        (gv).geom AS the_geom, 
        (gv).val
 FROM (SELECT id, 
              ST_Intersection(rast, the_geom) AS gv
       FROM srtm_tiled,
            cariboupoint_buffers_wgs
       WHERE ST_Intersects(rast, the_geom)
      ) foo;

Значения были затем обобщены с использованием следующего:

 CREATE TABLE result01 AS
 SELECT id, 
        sum(ST_Area(ST_Transform(the_geom, 32198)) * val) / 
        sum(ST_Area(ST_Transform(the_geom, 32198))) AS meanelev
 FROM caribou_srtm_inter
 GROUP BY id
 ORDER BY id;

На самом деле я не знаю достаточно PostGIS, чтобы объяснить это, но это звучит так, как будто вы пытались достичь. Учебник должен пролить свет на промежуточные этапы. Удачи :)

RK
источник
Спасибо @RK, я прочитал этот учебник. Я думаю, что мои расчеты более просты, но я все еще вычисляю этот базовый шаг!
djq
2

Что касается пункта 2 в первоначальном вопросе - в нескольких выпусках разработки Postgis 2.0 использовалась версия библиотеки GDAL, которая преобразует числа с плавающей точкой в ​​целые. Если в вашем растре есть значения с плавающей точкой, и вы использовали версию GDAL ниже 1.9.0 или предварительную версию PostGIS 2.0, которая неправильно вызывала GDALFPolygonize (), то вы можете столкнуться с этой ошибкой. Билеты в трекерах ошибок PostGIS и GDAL были зарегистрированы и закрыты. Эта ошибка была активна во время первоначального вопроса.

С точки зрения производительности вы обнаружите, что использование ST_Intersects(raster, geom)намного быстрее, чем использование ST_Intersects(geom, raster). Первая версия растеризует геометрию и выполняет пересечение растрового пространства. Вторая версия векторизует геометрию и делает пересечение векторного пространства, что может быть гораздо более дорогим процессом.

Пит Кларк
источник
0

Я также имеющие странные вопросы , используя ST_SummaryStatsс ST_Clip. Запросы к данным по-разному говорят мне, что минимальное значение моего растра было 32, а затем максимум 300, но ST_SummaryStatsвозвращалось значение -32700 для значений пикселей в моем целевом многоугольнике.

Я закончил взломом вокруг вопроса, как это:

WITH first AS (
   SELECT id, (ST_Intersection(geom, rast)).val
   FROM raster_table
   INNER JOIN vector_table ON ST_Intersects(rast, geom)
)
SELECT id, COUNT(val), SUM(val), AVG(val), stddev(val), MIN(val), MAX(val)
FROM first
GROUP BY id
jczaplew
источник