作品分享-用三十八行代码找到狭长面

内容简介

今天分享的这份作品，是作品征集发出后收到的第一份作品。这份作品，来自网名xiaofei的网友，现将作品整理后分享给大家。

作品介绍

我收到的作品，有处理后的文件，有源码，也有作者解答问题的思路。从形式上来说，非常的清晰，让人一看就能明白作者是怎么解答这个问题的。

以下为作品详细内容：

作品内容

要解答的问题

问题：如何找到并去除图斑的狭长面？

作者代码

library(rgeos)
library(rgdal)
library(magrittr)
library(spdep)
library(raster)

#自定义函数
splitTB<- function(x,width=60){
  coord<-x@polygons[[1]]@Polygons[[1]]@coords
  nb<-dnearneigh(coord,d1=1,d2=width,longlat = FALSE)
  al_line<-nb2lines(nb,coords = coordinates(coord))
  just<-gWithin(al_line,x,byid=T)
  fg_line<-al_line[c(just),]
  #线做buffer并于原始面做差，得到碎面
  line_buff<-gBuffer(fg_line,width = 0.000001)
  tb_sm<-gDifference(x,line_buff,byid=T) %>% disaggregate()
  tb_sm_id<-SpatialPolygonsDataFrame(tb_sm,data.frame(id=1:length(tb_sm)))
  #碎面质心做buff
  pts<-coordinates(tb_sm_id) %>% SpatialPoints()%>% SpatialPointsDataFrame(data.frame(ID=1:length(tb_sm_id)))
  pts_buff<-gBuffer(pts,width=width,byid=T)
  tb_ys<-spPolygons(coord) 
  pts_buff_dif<-(pts_buff - tb_ys)
  tc_id<- pts_buff_dif$ID 
  #结果输出
  tb_blpart<- tb_sm_id[-tc_id,] 
  tb_blpart$category <- 'BL'
  tb_tcpart<- (x-tb_blpart) 
  tb_tcpart$category <- 'TC'
  result<- rbind(tb_blpart,tb_tcpart)
}
#执行函数输出结果

shp<- readOGR("E:/Rdata/827/testnew.shp")
result<-lapply(1:nrow(shp), function(x){
  ls<- splitTB(shp[x,])
  }) %>% Reduce(rbind,.)
#写出结果
writeOGR(result,"E:/Rdata/827",layer='result827',driver = 'ESRI Shapefile')

解决思路：

通过R语言实现，采用大问题逐步拆解为小问题的思路；

1、如何找到宽度小于60m的部分？

宽度小于60m的部分，转化为找到相邻两点（空间位置上相邻）的间隔小于60m的部分，用R语言dnearneight找到相邻的两点，并用nb2lines输出两点连线。 作者在这里，将宽度小于60米的问题转换为点之间距离的计算问题，并且符合设置条件的点对连成了线。

2、从图斑中分割狭长部分?

第一想法是直接用线进行分割，但知识储备有限在R中没能实现，转而求其次，用线生成面宽度设置成0.000001，此容差对于arcgis来说几乎没有影响；然后再用原始图斑与线生成的缓冲面做差会将狭长部分与保留部分分开，狭长部分会被分割为多个碎面。 这一步中，作者采用了用线缓冲成面来切分另一个面的方式，来实现线分割面的功能。在主流的平台中，都没有直接用线分割面的功能，但这并不意味着不能用线来分割面，我们可以将面转成线，然后再重构面的方式来完成线分割面的功能，这样就不会有任何微小的缝隙。

3、如何将狭长部分生成的碎面识别出来？

通过buff面来识别：因为设置的规则是小于60m的狭长部分，所以就以第二步生成的图斑的质心为中心，做60m的缓冲面，如果60m的缓冲面超过图斑边界了那就说明宽度小于60m。 作者的这种方式很巧妙，通过缓冲面的分析，来找出狭长面，很棒！

以下内容就是在R语言中具体的实现方式了，我就不再评论了（主要我也不懂R）

3.1、第3步问题再细化：

① 如何知道缓冲面超过图斑边界了？ 问题转化为，用缓冲面跟图斑做差，如果缓冲面有剩余那就是超过边界了。 ② 如何根据缓冲面把碎图斑识别出来？ 通过赋id值的方式找出：先给第二步生成的图斑赋id，然后生成质心再赋id，然后对质心做缓冲面，这样缓冲面就带id了而且跟图斑id是对应关系，这样就知道哪些id的缓冲面是超边界了，然后根据id将碎图斑识别出来； ③ 对于有空洞的大图斑，生成的质心的缓冲面跟大图斑做差后也可能有剩余。 解决：用最外围的坐标点生成的图斑，这样就不会有空洞，做差后就不会有剩余。

4、输出狭长部分

用第二步生成的图斑将保留部分提取出来，然后再用原始图斑跟保留部分做差即可。对保留图斑和狭长部分分别赋上属性这样就可以在分别开来。 这里为什么不用碎图斑直接生成狭长面呢？ 因为在第二步的时候本想用线分割面，但是没实现，所以产生了0.00001的缝隙，无法直接融合为一块；

有待改进：

1、 后期有时间了研究一下怎么在R中用线直接分割面，用线分割面比做0.000001m的缓冲区更严谨。 2、 由于R容差的原因，导致生成狭长部分不够完美，有特别细的细缝，经实验比如s2在s1图斑内,s2边界正常，然后我用s1-(s1-s2)得到的图斑按常理说应该跟s2保持一致，但是现实确实多了细缝，我想这可能就是容差的问题，后期有时间了研究一下容差的设置；

最后输出结果图：

总结

总体来说，是一份很棒的作品不过作者可能再处理数据方面的经验稍微有点欠缺，导致做出的结果不是特别的完美。不过瑕不掩瑜，稍加改进，这份作品就会变得很完美。

由于篇幅限制，剩下的几份作品，我下次再分享。