Я пытаюсь вычислить растровую статистику (мин, макс, среднее) для каждого полигона в векторном слое, используя PostgreSQL / PostGIS.
В этом ответе GIS.SE описывается, как это сделать, путем вычисления пересечения между полигоном и растром, а затем вычисления средневзвешенного значения: https://gis.stackexchange.com/a/19858/12420.
Я использую следующий запрос (где dem
мой растр, topo_area_su_region
мой вектор и toid
уникальный идентификатор:
SELECT toid, Min((gv).val) As MinElevation, Max((gv).val) As MaxElevation, Sum(ST_Area((gv).geom) * (gv).val) / Sum(ST_Area((gv).geom)) as MeanElevation FROM (SELECT toid, ST_Intersection(rast, geom) AS gv FROM topo_area_su_region,dem WHERE ST_Intersects(rast, geom)) foo GROUP BY toid ORDER BY toid;
Это работает, но это слишком медленно. Мой векторный слой имеет 2489 тыс. Объектов, на обработку каждого из которых уходит около 90 мс - на обработку всего слоя уйдут дни . Скорость вычислений не кажется значительно улучшенной, если я вычисляю только min и max (что позволяет избежать вызовов ST_Area).
Если я выполняю аналогичные вычисления с использованием Python (GDAL, NumPy и PIL), я могу значительно сократить время, необходимое для обработки данных, если вместо векторизации растра (используя ST_Intersection) я растеризую вектор. Смотрите код здесь: https://gist.github.com/snorfalorpagus/7320167
Мне действительно не нужно взвешенное среднее - подход «если он касается, он в» достаточно хорош - и я достаточно уверен, что именно это замедляет ход событий.
Вопрос : есть ли способ заставить PostGIS вести себя так? т.е. возвращать значения всех ячеек из растра, к которому касается многоугольник, а не точное пересечение.
Я очень новичок в PostgreSQL / PostGIS, так что, возможно, есть что-то еще, что я не делаю правильно. Я использую PostgreSQL 9.3.1 и PostGIS 2.1 в Windows 7 (2,9 ГГц i7, 8 ГБ ОЗУ) и настроил конфигурацию базы данных, как указано здесь: http://postgis.net/workshops/postgis-intro/tuning.html
источник
Ответы:
Вы правы, использование
ST_Intersection
замедляет ваш запрос заметно.Вместо использования
ST_Intersection
лучше обрезать (ST_Clip
) ваш растр с помощью полигонов (ваших полей) и вывести результат как полигоны (ST_DumpAsPolygons
). Таким образом, каждая растровая ячейка будет преобразована в маленький многоугольник с разными значениями.Для получения минимального, максимального или среднего значения из дампов вы можете использовать те же операторы.
Этот запрос должен помочь:
В операторе
ST_Clip
вы определяете растр, растровую полосу (= 1), многоугольник и, если кадрирование должно быть ИСТИНА или ЛОЖЬ.Кроме того, вы можете использовать
avg((gv).val)
для расчета среднего значения.РЕДАКТИРОВАТЬ
Результат вашего подхода более точный, но медленный. Результаты объединения
ST_Clip
иST_DumpAsPolygons
игнорирования растровых ячеек, которые пересекаются с менее чем 50% (или 51%) их размера.Эти два снимка экрана с пересечения CORINE Land Use показывают разницу. Первая картинка с
ST_Intersection
, вторая сST_Clip
иST_DumpAsPolygons
.источник