微信公众号搜"智元新知"关注
微信扫一扫可直接关注哦!

如何计算 R 中不与其他多边形重叠的多边形面积?

如何解决如何计算 R 中不与其他多边形重叠的多边形面积?

我有一个 SpatialPoints 对象,其中包含多个兴趣点的坐标。 我还有几个形状文件(多边形),其中包含有关贫民窟存在的信息。每个 shapefile 中包含贫民窟信息的多边形可以重叠(它们提供的关于贫民窟存在的信息多少有些相同,但来自不同的来源)。

对于 SpatialPoints 对象中的每个点,我使用函数 spCircle 在每个点周围创建一个圆形多边形。接下来我需要做的是检查圆形多边形面积的多少百分比包含贫民窟。如果任何 shapefile 表明存在贫民窟,我会认为该地区存在贫民窟。

我创建了以下图像来帮助解释我的问题。圆圈表示围绕单个点的多边形。对于这一点,四个 shapefile 中的每一个都表明贫民窟存在于有些不同的区域(有时重叠,有时不重叠)。我希望能够找到红色区域(其中没有任何 shapefile 指示贫民窟的存在,然后计算有贫民窟的圆圈的百分比。

enter image description here

以下代码试图做到这一点:

# Create data with coordinates

lat = c(-22.868879628748203,-22.88511,-22.82166,-22.89692,-22.67945)
long = c(-43.237195000177564,-43.34278,-43.04717,-43.35168,-43.59667)

data_points = cbind.data.frame(lat,long)

coordinates(data_points) = c("lat","long")
proj4string(data_points) = CRS("+init=epsg:4326") 

# Transform projection of points to UTM

utmStr <- "+proj=utm +zone=%d +datum=NAD83 +units=m +no_defs +ellps=GRS80"
crs <- CRS(sprintf(utmStr,23))
data_points = spTransform(data_points,crs)


# Create a list with circular polygons around each point (radius = 2000 meters)

circular_grid = list()

for (i in 1:length(data_points)){
  
  spc = spCircle(radius = 2000,centerPoint = c(x=as.numeric(data_points@coords[i,1]),y=as.numeric(data_points@coords[i,2])),spID=i,spUnits = CRS("+proj=utm +zone=23 +datum=NAD83 +units=m +no_defs"))
  circular_grid[[i]] = spc
}


# For each circle,check the percentage that overlaps with several different shapefiles:

# I first use gUnion to merge all the shapefiles with info about slums together

allShapes = gUnion(shape1,shape2)
allShapes = gUnion(allShapes,shape3)
allShapes = gUnion(allShapes,shape4)
allShapes = gUnion(allShapes,shape5)
allShapes = gUnion(allShapes,shape6)
allShapes = as(allShapes,"SpatialpolygonsDataFrame")

allShapes = spTransform(allShapes,CRS("+proj=utm +zone=23 +datum=NAD83 +units=m +no_defs"))

# I am unable to reproduce the object "allShapes" (I do not kNow how),# but this is its information

# class       : SpatialpolygonsDataFrame 
# features    : 1 
# extent      : 633347.1,724692.1,-2547787,-2513212  (xmin,xmax,ymin,ymax)
# crs         : +proj=utm +zone=23 +datum=NAD83 +units=m +no_defs 
# variables   : 1
# names       : dummy 
# value       :     0

# Next,to get the intersection,I tried the following:

intersection_circle_shape = list()

for (i in 1:length(circular_grid)){
  
  circle = circular_grid[[i]][["spCircle"]]
  
  inter = intersect(circle,allShapes)
  
  intersection_circle_shape[[i]] = inter
  
}


# The list "intersection_circle_shape" is empty because the command
# "intersect" says that there is no intersection,but I kNow there is.

有什么想法吗?

版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。