Самый быстрый способ конвертировать большой растр в полилинию, используя R или Python?

14

У меня есть большой растровый файл (129600 на 64800 пикселей) с глобальными водоемами (значения 1 бит 0 и 1), и я пытаюсь извлечь береговые линии океана и внутренних вод.

Я пытался с ArcGIS и QGIS конвертировать из растра в полилинию, но это занимает много времени.

Кто-нибудь знает лучший / более быстрый способ (Python или R) или лучший инструмент для этой задачи?

Обновить

  • R: rasterToContour может быть быстрым и точным, но если у вас очень большой набор данных, такой как у меня (8 398 080 000 пикселей), вам нужен либо очень большой объем ОЗУ (более 16 ГБ), либо вы заставляете R делать больше обработки на жестком диске, и это также займет много времени.
  • Python / GDAL: gdal_poligonize создает полигоны вместо полилиний

Обновление 2

  • R rasterToContour: rasterToContour не дает желаемых результатов. По сравнению с ArcGIS (растр-полигон, за которым следует объект-линия), он не извлекает точный контур пикселя, как показано в примерах ниже.

Результат rasterToContour Результат rasterToContour

Результат ArcGIS Результат ArcGIS

ОБНОВЛЕНИЕ 3

Python / GDAL: я запустил gdal_polygonize из командной строки для ArcGIS на тестовом наборе данных, и результаты были предельно ясны:

  • гдал: 49 секунд
  • ArcGIS: 1,84 секунды
python r performance Общие Wevers
источник
Сделал это, см. Обновление 3.
Общие Wevers
Можете ли вы предоставить этот набор тестовых данных, чтобы мы могли видеть, являются ли предложенные альтернативы более быстрыми и / или дают ли требуемые результаты?
Керстен
Для такого огромного растра было бы лучше использовать C / C ++ с библиотекой gdal.
Родриго

Ответы:

8

Я работаю с R и использовал rasterToPolygonsиз rasterпакета в прошлом, но теперь я предпочитаю gdal_polygonizeRДжона Баумгартнера. Он основан на gdal_polygonize.pyи намного быстрее. Джон Баумгартнер опубликовал код и привел пример для использования в своем блоге .

Если вы знакомы с Python, вы можете использовать gdal_polygonize.pyнапрямую, конечно.

радужная оболочка
источник
1
Я дам это, я пытаюсь. В прошлый раз я использовал gdal_polygonize.py ArcGIS был еще быстрее.
Общие Wevers
Я не ожидал, что ArcGis может быть быстрее, чем GDAL. @ Generic Militzer
Iris,
Ах, подождите, это создаст полигоны, но мне нужны полилинии.
Общие Wevers
Если вы поместите свои данные в базу геоданных, это будет довольно быстро. Но все еще недостаточно быстро. Вот почему я ищу альтернативы.
Общий Wevers
2
Это не обязательно проблема получения полигонов, вы всегда можете преобразовать их в полилинии (хотя при таком количестве это, конечно, может занять некоторое время).
Мартин
7

Для потомков, у меня был успех с stars::пакетом Rдля выполнения этого типа операции быстро.

library(raster)
library(stars)
library(sf)
library(magrittr)

f <- system.file("external/test.grd", package="raster")
r <- raster(f)
r[r[] < 750] <- 0
r[r[] >= 750] <- 1

x <- st_as_stars(r) %>% 
  st_as_sf() %>% # this is the raster to polygons part
  st_cast("MULTILINESTRING") # cast the polygons to polylines

plot(x)

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

plot(r)
plot(x, add = TRUE)

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

mikoontz
источник
5

Попробуйте rasterToContourиз растрового пакета.

f <- system.file("external/test.grd", package="raster")
r <- raster(f)
r[r[] < 750] <- 0
r[r[] >= 750] <- 1

x <- rasterToContour(r)
class(x)
> [1] "SpatialLinesDataFrame"
> attr(,"package")
> [1] "sp"

plot(r)
plot(x, add=TRUE)

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

Затем вы можете легко записать файлы в локальную папку, например, как «ESRI Shapefile» (.shp), используя приведенный ниже код. Посмотрите ogrDriversиз rgdal , чтобы выяснить , какие драйверы ваша система совместима с.

library(rgdal)
writeOGR(x, dsn = getwd(), layer = "coastlines", driver = "ESRI Shapefile")
fdetsch
источник
Я постараюсь держать пальцы скрещенными, это не убьет мою оперативную память. Несмотря на то, что у меня есть 16 ГБ, чего, надеюсь, достаточно, R иногда не так эффективен с большими растровыми файлами. Но посмотрим.
Общие Wevers
Конверсия работала как-то, но я не смог проверить подробно. Поскольку я обычно больше занимаюсь обработкой растровых данных, вы можете сказать мне, как я могу перенести SpatialLineDataFrame в шейп-файл или что-то сопоставимое. Я гуглил и все еще борюсь, так как не знаю названия слоя (OGRwrite).
Общий Wevers
Ха-ха, я определенно понимаю твою точку зрения. Смотрите обновление выше.
fdetsch
2
Еще один совет: попробуйте установить для параметра maxpixels rasterToContourболее высокое значение, например 1e + 9. В итоге вы получите больше деталей. Настройка по умолчанию создает довольно обобщенные контурные линии.
fdetsch
1
Если вы не хотите, чтобы resampleваши данные имели более грубое пространственное разрешение, единственное решение, которое я могу себе представить, это разделить ваши данные на несколько плиток (например, 16 вспомогательных растров), а затем выполнить rasterToContourитерацию для каждой плитки отдельно и и, наконец, mergeполучившиеся шейп-файлы в один огромный шейп-файл. Если вас это интересует, пакет нашей рабочей группы Rsenal предлагает функцию, которая вызываетsplitRaster создание нескольких вспомогательных растров из одного огромного растра.
fdetsch
2

Хотя я большой поклонник GDAL, инструмент полигонизации был слишком медленным для моих приложений.

Быстрая альтернатива - gdal_trace_outlineот сценариев Dans GDAL, которые также имеют больше вариантов относительно толерантности, пончиков и т. Д.

Подобно gdal_polygonizeэтому, также создаются полигоны, которые потом нужно будет конвертировать ogr2ogr -nlt MULTILINESTRING.

Недостатком является то, что вам нужно скомпилировать его самостоятельно, если вы не используете Linux или Mac OsX System.

Керстен
источник
К сожалению, это не удалось с сообщением об ошибке: «Ошибка сегментации (ядро сброшено)». Я предполагаю, что мой файл слишком большой или более точный, он будет производить слишком много маленьких полигонов.
Общий Wevers