Создание большого количества случайных точек в двоичном растре?

Я хочу создать точечный векторный набор данных из 10000 точек (или больше) внутри двоичного растра, где точки должны быть ограничены областями, где значение растра равно 1.

Я попробовал следующие шаги.

1. Многоугольный растр
2. QGIS: Вектор -> Инструменты исследования -> Случайные точки

Это прекрасно работает до 2000 очков, но все, что выше, просто вызывает сбой QGIS.

Есть ли способ создать набор векторных данных с большим количеством точечных объектов, ограниченных двоичным растром (или его полигональной версией)?

Следующие инструменты находятся в моем распоряжении, ранжируются от самых неблагоприятных: QGIS, Python, R, ArcGIS

Это то, к чему я стремлюсь, только с 10-кратными точечными функциями.

Вот способ в R:

Создайте тестовый растр, 20х30 ячеек, сделайте 1/10 из ячеек, установленных на 1, построите график:

> require(raster)
> m = raster(nrow=20, ncol=30)
> m[] = as.numeric(runif(20*30)>.9)
> plot(m)

Для существующего растра в файле, например, geoTIFF, вы можете просто сделать:

> m = raster("mydata.tif")

Теперь возьмите матрицу координат xy ячеек 1, нарисуйте эти точки, и мы увидим, что у нас есть центры клеток:

> ones = xyFromCell(m,1:prod(dim(m)))[getValues(m)==1,]
> head(ones)
       x    y
[1,] -42 85.5
[2,] 102 85.5
[3,] 162 85.5
[4,]  42 76.5
[5,] -54 67.5
[6,]  30 67.5
> points(ones[,1],ones[,2])

Шаг 1. Создайте 1000 (xo, yo) пар, которые центрированы на 0 в коробке размером с одну ячейку. Обратите внимание на использование, resчтобы получить размер ячейки:

> pts = data.frame(xo=runif(1000,-.5,.5)*res(m)[1], yo=runif(1000,-.5,.5)*res(m)[2])

Шаг 2. Определите, в какую ячейку входит каждая из вышеперечисленных точек, путем случайной выборки 1000 значений от 1 до числа 1:

> pts$cell = sample(nrow(ones), 1000, replace=TRUE)

Наконец, вычислите координату, добавив центр ячейки к смещению. Участок для проверки:

> pts$x = ones[pts$cell,1]+pts$xo
> pts$y = ones[pts$cell,2]+pts$yo
> plot(m)
> points(pts$x, pts$y)

Вот 10000 точек (замените 1000 выше на 10000), построенных с помощью pch=".":

Практически мгновенно для 10000 точек на растре 200x300 с половиной точек как единица. Думаю, со временем они будут линейно увеличиваться с тем, сколько в растре.

Чтобы сохранить как шейп-файл, преобразуйте его в SpatialPointsобъект, задайте ему правильную ссылку на систему координат (так же, как и в растре) и сохраните:

> coordinates(pts)=~x+y
> proj4string(pts)=CRS("+init=epsg:4326") # WGS84 lat-long here
> shapefile(pts,"/tmp/pts.shp")

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

Всякий раз, когда я работаю с большими наборами данных, мне нравится запускать инструменты / команды вне QGIS, например, из автономного скрипта или из оболочки OSGeo4W . Не столько потому, что QGIS дает сбой (даже если он говорит «Не отвечает», возможно, он все еще обрабатывает данные, которые вы можете проверить из диспетчера задач ), а потому, что для обработки данных доступно больше ресурсов ЦП, таких как ОЗУ. Сам QGIS потребляет изрядную часть памяти для запуска.

В любом случае, чтобы запустить инструмент вне QGIS ( вам необходимо установить QGIS через установщик OSGeo4W ), выполните первые 2 шага, как описано @gcarrillo в этом посте: Проблема с импортом qgis.core при написании автономного скрипта PyQGIS (Я предлагаю скачать и использовать его .bat файл).

Как только ПУТИ установлены, введите pythonв командной строке. Для удобства скопируйте следующий код в текстовый редактор, такой как «Блокнот», отредактируйте параметры, такие как путь к вашему шейп-файлу и т. Д., А затем вставьте все это в командную строку, щелкнув правой кнопкой мыши> Вставить :

import os, sys
from qgis.core import *
from qgis.gui import *
from PyQt4.QtGui import *

from os.path import expanduser
home = expanduser("~")

QgsApplication( [], False, home + "/AppData/Local/Temp" )

QgsApplication.setPrefixPath("C://OSGeo4W64//apps//qgis", True)
QgsApplication.initQgis()
app = QApplication([])

sys.path.append(home + '/.qgis2/python/plugins')
from processing.core.Processing import Processing
Processing.initialize()
from processing.tools import *

shape = home + "/Desktop/Polygon.shp"
result = home + "/Desktop/Point.shp"
general.runalg("qgis:randompointsinlayerbounds", shape, 10000, 0, result)

Используя скрипт, я запустил инструмент Случайные точки в границах слоев для довольно большого шейп-файла, и для генерации 10 тыс. Точек потребовалось менее 20 секунд. Запуск его внутри QGIS занял почти 2 минуты, так что, по крайней мере, для меня есть существенная разница.

Вы также можете использовать GRASS GIS непосредственно для этой работы - Стратифицированная случайная выборка: случайная выборка из векторной карты с пространственными ограничениями :

https://grass.osgeo.org/grass72/manuals/v.random.html#stratified-random-sampling:-random-sampling-from-vector-map-with-spatial-constraints

Кроме того, в команде реализована случайная выборка из векторной карты по атрибуту и ​​несколько других методов.

Примечание. Версия v.random, представленная в QGIS посредством обработки, отражает не полную функциональность, а просто упрощенное представление.