这周我一直在玩terra中的buffer函数。我应用了一个负缓冲区来获得我研究区域内的一个区域。我最初的输入多边形创建了多个奇怪的结果，并多次崩溃R。我发现了以下解决方案，但我怀疑它们是变通方案，可能有更好的方法来解决它们。

• 问题1：聚合多边形会在原始多边形之外创建一个奇数多边形
• 解决方案：在原始文件上运行disagg
• 问题2：带有孔（或内环）的多边形在孔/环的外部创建了一个缓冲区，该缓冲区可以超出外部多边形
• 解决方案：在原始文件上运行fillHoles
• 问题3：R经常崩溃（尽管不规则）
• 解决方法：使用非投影坐标系代替输入文件的原始投影坐标系

我上周更新了terra，在Windows 11上运行的是1.7-55版本。我有一个第13代英特尔（R）酷睿（TM）i7- 1355 U 1.70 GHz和16 GB RAM。

library(terra)

issue 1 : aggregated polygons create an odd polygon that is outside the original polygons

#three-polygon spatVector which generates odd inset buffer
p1<-vect("POLYGON ((-66.940374 45.008097, -66.971762 45, -66.986546 45.008447, -66.940267 45.026992, -66.940374 45.008097))")                                                                                                                                                                                                              
p2<-vect("POLYGON ((-66.828 45.032592, -66.829 45.035, -66.824 45.031, -66.8301 45.0275, -66.834 45.027, -66.828 45.032592))")                                                                                                                                                                                 
p3<-vect("POLYGON ((-67 46, -66.2 46, -66.2 45.65, -66.37 45.45, -66.37 45.105, -66.459 45.05, -66.46 45.17,
-66.78 45.04, -66.76 45.12, -66.89 45.04, -66.916 45.11, -66.89 45.1, -66.82 45.127, -66.96 45.19,
-67 45.15, -67 46))")

wkt.p<-rbind(p1, p2, p3)
crs(wkt.p)<-"+proj=longlat +datum=WGS84"

hol_buff<-buffer(aggregate(wkt.p), -2000)
plot(wkt.p, ylim=c(45, 46))
plot(hol_buff, add=T, col="green")

#dissaggreated polygons do not create odd inset
hol_buff<-buffer(wkt.p, -2000)
plot(hol_buff, add=T, col="red")

#three-polygon spatVector which does NOT generates odd inset buffer (upper limit of p3 changed from 46 to 45.8)
p1<-vect("POLYGON ((-66.940374 45.008097, -66.971762 45, -66.986546 45.008447, -66.940267 45.026992, -66.940374 45.008097))")                                                                                                                                                                                                              
p2<-vect("POLYGON ((-66.828 45.032592, -66.829 45.035, -66.824 45.031, -66.8301 45.0275, -66.834 45.027, -66.828 45.032592))")                                                                                                                                                                                 
p3<-vect("POLYGON ((-67 45.8, -66.2 45.8, -66.2 45.65, -66.37 45.45, -66.37 45.105, -66.459 45.05, -66.46 45.17,
-66.78 45.04, -66.76 45.12, -66.89 45.04, -66.916 45.11, -66.89 45.1, -66.82 45.127, -66.96 45.19,
-67 45.15, -67 45.8))")

wkt.p<-rbind(p1, p2, p3)
crs(wkt.p)<-"+proj=longlat +datum=WGS84"

hol_buff<-buffer(aggregate(wkt.p), -2000)
plot(wkt.p, ylim=c(45, 46))
plot(hol_buff, add=T, col="blue")

字符串
x1c 0d1x的数据

#issue no 2 : holes within polygons leads to negative buffer that are outside the polygon boundary

par(mfrow=c(1,2))
#polygon with inner ring
wkt.y<-"POLYGON ((-66.902 45.013, -66.901 45.01, -66.906 45.0114,  -66.9045 45.0155, -66.902 45.0155, -66.902 45.013),
(-66.903371 45.012734, -66.903372 45.012914, -66.903626 45.012913, -66.903624 45.012733, -66.903371 45.012734))"
y<-vect(wkt.y,crs="+proj=longlat +datum=WGS84")
plot(y, main="terra output")
plot(buffer(y, -40), add=T, col="red")

#polygon without inner ring
wkt.w<-"POLYGON ((-66.902 45.013, -66.901 45.01, -66.906 45.0114,  -66.9045 45.0155, -66.902 45.0155, -66.902 45.013))"
w<-vect(wkt.w,crs="+proj=longlat +datum=WGS84")
plot(y, main="desired output")
plot(buffer(w, -40), add=T, col="red")#negative buffer using the simple polygon that shows the desired buffered polygon
plot(y, add=T)#hypothetical input
#removing the inner ring can be done by using fillHoles
w2<-fillHoles(y)

#polygon with inner ring at edge
par(mfrow=c(1,2))
wkt.y<-"POLYGON ((-66.902 45.012,  -66.906 45.0114,  -66.906 45.0155, -66.902 45.0155, -66.902 45.012),
(-66.903371 45.012734, -66.903372 45.012914, -66.903626 45.012913, -66.903624 45.012733, -66.903371 45.012734))"
y<-vect(wkt.y,crs="+proj=longlat +datum=WGS84")
#plot(buffer(y, -200), col="red", )
plot(y, xlim=c(-66.907, -66.901), ylim=c(45.01, 45.017), main="terra output")
plot(buffer(y, -200), add=T, col="red", border="forestgreen", lwd=2)
plot(y, add=T)#overlay original polygon

#polygon without inner ring
wkt.w<-"POLYGON ((-66.902 45.012,  -66.906 45.0114,  -66.906 45.0155, -66.902 45.0155, -66.902 45.012))"
w<-vect(wkt.w,crs="+proj=longlat +datum=WGS84")
plot(y, main="desired output")
plot(buffer(w, -200), add=T, col="red")#negative buffer using the simple polygon that shows the desired buffered polygon
plot(y, add=T)#hypothetical input

型