#r #raster #spatial
Вопрос:
Я несколько дней боролся за создание растра, отображающего расстояние до дорог. У меня есть файл формы пространственных линий дорог в моем регионе исследования, и я хочу создать объект «расстояние до дороги» с тем же разрешением, что и другие мои ковариаты, 30-метровые пиксели. За последние несколько дней я перепробовал несколько вещей, но совсем недавно я превратил свой растровый шаблон в пространственные точки и использовал gDistance , но, похоже, это продолжается и продолжается без конца.
Мой ландшафт несколько велик, но не огромен (~25000 км2), другой пространственный анализ не занимает так много времени, даже с большими стопками из 100 растров. На данный момент я несколько растерян, перепробовав большинство предложений, которые я нашел в Интернете.
Комментарии:
1. Можете ли вы предоставить данные строки или что-то в аналогичном формате? Можно ли проектировать в плоские координаты, а затем работать с пакетом spatstat?
2. Конечно, я могу предоставить данные о линии. Я не использовал пакет spatstat, но, похоже, все должно быть в порядке. Однако я не уверен, что лучший способ поделиться данными здесь, в переполнении стека.
3. В приведенном ниже ответе я просто использовал примеры данных из
sfпакета в R. Чтение данных в строке и преобразование в формат spatstat иногда происходит мучительно медленно, если у вас много строк. Возможно, прохождениеas.linnet.SpatialLinesmaptoolsможет ускорить процесс, и, возможно, вы могли бы упростить дорожную сеть, если она будет очень подробной.
Ответ №1:
Чтение данных шейп-файла с sf помощью (обратите внимание, что пользовательская система координат CRS, используемая для проецирования в плоское пространство, была найдена с использованием https://projectionwizard.org):
library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
dat_path <- system.file("shape/olinda1.shp", package="sf")
dat <- st_read(dat_path)
#> Reading layer `olinda1' from data source `/usr/lib/R/site-library/sf/shape/olinda1.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 470 features and 6 fields
#> Geometry type: POLYGON
#> Dimension: XY
#> Bounding box: xmin: -34.91692 ymin: -8.044467 xmax: -34.82779 ymax: -7.954672
#> Geodetic CRS: GRS 1980(IUGG, 1980)
poly <- st_geometry(dat)
new_crs <- st_crs(" proj=tcea lon_0=-34.87 datum=WGS84 units=m no_defs")
flat_poly <- st_transform(poly, new_crs)
flat_lines <- st_boundary(flat_poly)
plot(flat_lines)
Преобразование в формат spatatstat (обратите внимание, что spatstat не выполняет автоматического
преобразования единиц измерения или именования, поэтому эта часть кода предназначена только для удобного отображения
. Кроме того, координаты смещены, поэтому начало координат находится в нижней части
для более качественного вывода на печать, но если вы когда-либо конвертируете обратно в другие форматы, это
следует опустить, чтобы избежать неправильного расположения области на карте.)
library(spatstat)
#> Loading required package: spatstat.data
#> Loading required package: spatstat.geom
#> spatstat.geom 2.2-2.002
#> Loading required package: spatstat.core
#> Loading required package: nlme
#> Loading required package: rpart
#> spatstat.core 2.3-0.003
#> Loading required package: spatstat.linnet
#> spatstat.linnet 2.3-0
#>
#> spatstat 2.2-0 (nickname: 'That's not important right now')
#> For an introduction to spatstat, type 'beginner'
flat_lines_psp <- as.psp(flat_lines)
flat_lines_psp
#> planar line segment pattern: 12235 line segments
#> window: rectangle = [-5172.766, 4653.619] x [-895506, -885510.1] units
L <- shift(flat_lines_psp, origin = "bottomleft")
unitname(L) <- "m"
L
#> planar line segment pattern: 12235 line segments
#> window: rectangle = [0, 9826.385] x [0, 9995.897] m
Растровое изображение по умолчанию с разрешением 128х128 пикселей:
system.time(D128pix <- distmap(L))
#> user system elapsed
#> 0.527 0.000 0.527
D128pix
#> real-valued pixel image
#> 128 x 128 pixel array (ny, nx)
#> enclosing rectangle: [0, 9826.4] x [0, 9995.9] m
plot(D128pix, log = TRUE)
Вы можете управлять разрешением карты расстояний с помощью любого dimyx
(количество результирующих пикселей) или eps (разрешение/размер пикселей):
system.time(D512pix <- distmap(L, dimyx = 512))
#> user system elapsed
#> 8.083 0.020 8.104
plot(D512pix, log = TRUE)
system.time(D20m <- distmap(L, eps = 20))
#> user system elapsed
#> 7.704 0.000 7.704
plot(D20m, log = TRUE)
Я думаю, что код растет более или менее линейно по количеству пикселей и
по количеству сегментов линий, поэтому, если у вас много небольших сегментов линий, это
может занять очень много времени. По крайней мере, в этом примере уменьшение количества
строк примерно в 2 раза сокращает время вычисления в 2 раза:
Lsmall <- L[owin(c(2500,7500),c(2500,7500))]
Lsmall
#> planar line segment pattern: 5988 line segments
#> window: rectangle = [2500, 7500] x [2500, 7500] m
system.time(Dsmall512pix <- distmap(Lsmall, dimyx = 512))
#> user system elapsed
#> 4.179 0.004 4.184
plot(Dsmall512pix, log = TRUE)
