쉐이프 파일에서 커버 백분율 계산하기

Question

나는 보통 쉐이프 파일과 함께 작동하지 않으므로 여기서 약간 잃어 버렸습니다. 두 개의 shapefile 각각에 여러 개의 객체가 있습니다. 첫 번째는 32 개의 다각형 세트입니다 (각각 하나의 플롯입니다). 두 번째 shapefile은 각 플롯 내에서 서로 다른 크기의 초목 군집을 나타내는> 10,000 개의 객체를 포함합니다. 나는 알아 내려고 애 쓰고있다.

1) 각 사이트 내 총 식생 커버의 비율을 어떻게 계산합니까?

2) 각 플롯에서 식물 덮개의 면적 중 5 % 미만이 몇 퍼센트입니까?

내 그림은 ArcGIS에서 단일 플롯의 데이터입니다.

Answer 1

다음 코드는 당신이 원하는 것을 할 것입니다, 나는 생각한다.

NB : 이는 아래 설명 된대로 shapefile 폴리곤에 저장된 영역 정보를 사용합니다. 그것은 가 아니며 식물의 shapefile 자료 섹션에있는 Area 열을 사용하십시오. 대부분의 경우 Area은 shapefile에 저장된 영역과 동일하지만 경우에 따라 Area이 훨씬 더 큽니다. 귀하의 Area 데이터가 어디서 왔는지 모르기 때문에 shapefile 다각형과 함께 저장된 정보를 사용하는 것이 더 안전합니다.

library(rgdal) 
library(ggplot2) 

setwd("<directory containing all your shapefiles>") 
plt.map <- readOGR(dsn=".",layer="plots") 
veg.map <- readOGR(dsn=".",layer="veg_in_plots") 
# associate LocCode with polygon IDs 
plt.data <- cbind(id=rownames(plt.map@data), LocCode=plt.map@data$LocCode) 
veg.data <- cbind(id=rownames(veg.map@data), LocCode=veg.map@data$LocCode) 
# function to extract area from polygon data 
get.area <- function(polygon) { 
    row <- data.frame(id=polygon@ID, area=polygon@area, stringsAsFactors=F) 
    return(row) 
} 
# area of each plot polygon 
plt.areas <- do.call(rbind,lapply(plt.map@polygons, get.area)) 
plt.data <- merge(plt.data,plt.areas, by="id") # append area column to plt.data 
# area of each vegetation polygon 
veg.areas <- do.call(rbind,lapply(veg.map@polygons, get.area)) 
veg.data <- merge(veg.data,veg.areas, by="id") # append area column to veg.data 
# total area of vegetation polygons by LocCode 
veg.smry <- aggregate(area~LocCode,data=veg.data,sum) 
smry  <- merge(plt.data,veg.smry,by="LocCode") 
smry$coverage <- with(smry,100*area.y/area.x) # coverage percentage 
# total area for vegetation object with A < 5 msq 
veg.lt5 <- aggregate(area~LocCode,data=veg.data[veg.data$area<5,],sum) 
smry  <- merge(smry, veg.lt5, by="LocCode") 
# fraction of covered area coming from veg. obj. with A < 5 msq 
smry$pct.lt5 <- with(smry, 100*area/area.y) 

이를 생성합니다 :

# head(smry) 
# LocCode id area.x area.y coverage  area pct.lt5 
# 1  1 3 1165.916 259.2306 22.23408 60.98971 23.52720 
# 2  10 11 1242.770 366.3222 29.47626 88.21827 24.08216 
# 3  11 12 1181.366 213.2105 18.04779 129.21612 60.60496 
# 4  12 13 1265.352 577.6037 45.64767 236.83946 41.00380 
# 5  13 14 1230.662 226.2686 18.38593 48.09509 21.25575 
# 6  14 15 1274.538 252.0577 19.77640 46.94874 18.62619 

설명 :

모양 파일이 rgdal 패키지에 readOGR(...)를 사용하여 R로 가져올 수 있습니다. 다각형 모양 파일을 가져올 때 결과는 "SpatialPolygonDataFrame"개체입니다. 이러한 객체는 기본적으로 두 개의 섹션으로 구성됩니다. 각 폴리곤을 그릴 때 필요한 좌표가있는 다각형 섹션과 각 다각형 (즉, 다각형 당 한 행)에 대한 데이터가있는 데이터 섹션이 있습니다. Shape 파일이 map 예컨대로 가져온 경우

map <- readOGR(dsn=".",layer="myShapeFile") 

는 폴리곤 데이터 섹션 및 map@polygonmap@data로 액세스 할 수있다. 폴리곤 영역은 폴리곤 영역에 저장됩니다. 영역을 얻으려면 폴리곤에서 영역과 폴리곤 ID를 추출하는 함수 get.area(...)을 정의합니다. 그런 다음 우리는 lapply(...)를 사용하여 모든 다각형에 대한 그 함수를 호출하고, 모든 반환 된 값이 함께 행 방향 rbind(...)을 사용하여 바인딩 :

plt.areas <- do.call(rbind,lapply(plt.map@polygons, get.area)) 
veg.areas <- do.call(rbind,lapply(veg.map@polygons, get.area)) 

이제 우리는 플롯 다각형과 식물의 영역을 연결해야합니다. 이 작업은 각 모양 파일의 데이터 섹션에있는 열 LocCode을 통해 수행됩니다.

plt.data <- cbind(id=rownames(plt.map@data), LocCode=plt.map@data$LocCode) 
veg.data <- cbind(id=rownames(veg.map@data), LocCode=veg.map@data$LocCode) 

그런 다음 우리는 다각형의 ID를 기반으로 지역의 열을 추가합니다 :

plt.data <- merge(plt.data,plt.areas, by="id") # append area column to plt.data 
veg.data <- merge(veg.data,veg.areas, by="id") # append area column to veg.data 

그런 다음 우리가 LocCode에 의해 식물의 영역을 합산 할 필요가 플롯과 식물 영역 모두 LocCode와 그래서 우리 첫째 준 다각형 ID의의 :

veg.smry <- aggregate(area~LocCode,data=veg.data,sum) 

그리고 마지막으로 플롯 다각형 영역으로이 병합 :

smry  <- merge(plt.data,veg.smry,by="LocCode") 

smry 데이터 프레임에서 area.x은 플롯의 영역이고 area.y은 해당 플롯에서 초목으로 덮인 총 영역입니다. 두 셰이프 파일 모두 투영법은 다음과 같습니다.

+proj=utm +zone=13 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0 

단위는 미터이고 영역은 msq입니다. 지역에서 < 5 MSQ 오는 방법 식물의 대부분을 결정하기 위해, 우리는 지역 < 5 식생 지역을 총과 smry와 그 결과를 병합 :

veg.lt5 <- aggregate(area~LocCode,data=veg.data[veg.data$area<5,],sum) 
smry  <- merge(smry, veg.lt5, by="LocCode") 

마지막으로, 데이터를 우리가지도를 렌더링하는 간단합니다 있습니다 각 플롯 영역 :

cols <- c("id","LocCode") 
plt.df <- fortify(plt.map) 
plt.df <- merge(plt.df,plt.data[cols],by="id") 
veg.df <- fortify(veg.map) 
veg.df <- merge(veg.df,veg.data[cols],by="id") 
ggp <- ggplot(plt.df, aes(x=long, y=lat, group=group)) 
ggp <- ggp + geom_path() 
ggp <- ggp + geom_polygon(data=veg.df, fill="green") 
ggp <- ggp + facet_wrap(~LocCode,scales="free") 
ggp <- ggp + theme(axis.text=element_blank()) 
ggp <- ggp + labs(x="",y="") 
ggp <- ggp + coord_fixed() 
ggp