问题描述
我正在尝试使用 R 中的 raster
包从光栅对象中提取轮廓线。
rasterToContour
似乎运行良好,绘图也很好,但经过调查,轮廓线似乎被分解成不规则的部分。来自 ?rasterToContour
library(raster)
f <- system.file("external/test.Grd",package="raster")
r <- raster(f)
x <- rasterToContour(r)
class(x)
plot(r)
plot(x,add=TRUE)
我正在尝试提取栅格中示例站点的轮廓线。因此,我们随机选择一个站点,提取其高程,然后再次运行 rasterToContour()
,指定等高线 level
的高程。
# our sample site - a random cell chosen on the raster
xyFromCell(r,5000) %>%
SpatialPoints(proj4string = crs(r)) %>%
{. ->> site_sp} %>%
st_as_sf %>%
{. ->> site_sf}
# find elevation of sample site,and extract contour lines
extract(r,site_sf) %>%
{. ->> site_elevation}
# extract contour lines
r %>%
rasterToContour(levels = site_elevation) %>%
{. ->> contours_sp} %>%
st_as_sf %>%
{. ->> contours_sf}
# plot the site and new contour lines (approx elevation 326)
plot(r)
plot(contours_sf,add = TRUE)
plot(site_sf,add = TRUE)
# plot the contour lines and sample site - using sf and ggplot
ggplot()+
geom_sf(data = contours_sf)+
geom_sf(data = site_sf,color = 'red')
然后我们使用 st_intersects
找到与站点相交的线(缓冲区宽度为 100 以确保它与线接触)。但是,这将返回所有等高线。
contours_sf %>%
filter(
st_intersects(.,site_sf %>% st_buffer(100),sparse = FALSE)[1,]
) %>%
ggplot()+
geom_sf()
我假设所有行都被返回,因为它们看起来是一个 MULTILInesTRING
sf
对象。
contours_sf
# Simple feature collection with 1 feature and 1 field
# geometry type: MULTILInesTRING
# dimension: XY
# bBox: xmin: 178923.1 ymin: 329720 xmax: 181460 ymax: 333412.3
# CRS: +proj=sterea +lat_0=52.1561605555556 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +datum=wgs84 +units=m +no_defs
# level geometry
# C_1 326.849822998047 MULTILInesTRING ((179619.3 ...
因此,我使用 contours_sf
将 MULTILInesTRING
ngeo::st_segments()
拆分为单独的行(我找不到任何 sf
方法来执行此操作,但我愿意使用替代方法,特别是如果这是问题所在)。
没想到这会返回 394 个特征;从看图来看,我预计大约有 15 行。
contours_sf %>%
nngeo::st_segments()
# Simple feature collection with 394 features and 1 field
# geometry type: LInesTRING
# dimension: XY
# bBox: xmin: 178923.1 ymin: 329720 xmax: 181460 ymax: 333412.3
# CRS: +proj=sterea +lat_0=52.1561605555556 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +datum=wgs84 +units=m +no_defs
# First 10 features:
# level result
# 1 326.849822998047 LInesTRING (179619.3 329739...
# 2 326.849822998047 LInesTRING (179580 329720.4...
# 3 326.849822998047 LInesTRING (179540 329720,...
# 4 326.849822998047 LInesTRING (179500 329735.8...
# 5 326.849822998047 LInesTRING (179495.3 329740...
# 6 326.849822998047 LInesTRING (179460 329764,...
# 7 326.849822998047 LInesTRING (179442.6 329780...
# 8 326.849822998047 LInesTRING (179420 329810,...
# 9 326.849822998047 LInesTRING (179410.2 329820...
# 10 326.849822998047 LInesTRING (179380 329847.3...
然后,当我们过滤以仅保留与站点相交的线(缓冲区宽度为 100)时,仅返回预期轮廓线的一小部分(线的红色部分,我假设反映了 100 缓冲区宽度)。
contours_sf %>%
nngeo::st_segments() %>%
filter(
# this Syntax used as recommended by this answer https://stackoverflow.com/a/57025700/13478749
st_intersects(.,sparse = FALSE)
) %>%
ggplot()+
geom_sf(colour = 'red',size = 3)+
geom_sf(data = contours_sf)+
geom_sf(data = site_sf,colour = 'cyan')+
geom_sf(data = site_sf %>% st_buffer(100),colour = 'cyan',fill = NA)
任何人都对以下几点有想法:
- 解释为什么等高线“断了”
- 提供一种将碎片“连接”在一起的有效方法
- 替代
nngeo::st_segments()
,如果这实际上是 394 行而不是 ~15 行的来源
解决方法
将 MULTILINESTRING 转换为 LINESTRING 似乎可以满足您的需求:
contours_sf %>% st_cast("LINESTRING") %>%
filter(st_intersects(.,st_buffer(site_sf,100),sparse=FALSE)[,1]) %>%
ggplot()+
geom_sf(data = contours_sf)+
geom_sf(data = site_sf,color = 'red') +
geom_sf(color = 'pink')
,
如果你从分解线条开始,也许效果会更好
dangerouslySetInnerHTML
或者用library(raster)
f <- system.file("external/test.grd",package="raster")
r <- raster(f)
x <- rasterToContour(r)
x <- disaggregate(x)
terra
你可以继续这样
library(terra)
r <- rast(f)
x <- as.contour(r)
x
# class : SpatVector
# geometry : lines
# dimensions : 8,1 (geometries,attributes)
x <- disaggregate(x)
x
# class : SpatVector
# geometry : lines
# dimensions : 56,attributes)
或者像这样
y <- st_as_sf(x)