问题描述
我正在尝试计算美国每个县的土地覆盖重新划分。 我已经使用 FedData 包(devtools 版本)获得了 Apache 县的 NLCD,我正在使用人口普查局的县 shapefile。
问题是我得到的区域比官方的和我的 shapefile 中指示的区域大得多,即 51,000km^2 而不是官方的 29,0000km^2。一定有一些我不了解光栅对象的地方,但经过数小时的网络搜索后,我感到非常困惑,感谢任何帮助。
下面介绍使用的代码和计算方法。县数据可以在这里下载: https://www2.census.gov/geo/tiger/TIGER2016/COUNTY/
以下代码假设县 shapefile 已保存并解压缩。
- 获取和读取数据
#devtools::install_github("ropensci/FedData")
library(FedData)
library(rgdal)
library(dplyr)
#Get Apache polygone
counties<- readOGR('tl_2016_us_county/tl_2016_us_county.shp')
apache <- subset(counties,counties$GEOID=="04001")
# Get NCLD data
nlcd_data <- get_nlcd(apache,year = 2011,label = "Apache",force.redo = TRUE)
nlcd_data #inspect the object,can see that number of cells is around 57 million
- 然后我提取了栅格的值并将它们放入频率表中。从那里我计算结果面积。由于 NLCD 数据是 30m 分辨率,我将每个类别的单元格数乘以 900,然后将结果除以 100 万,得到以平方公里为单位的面积。
计算的总面积过大。
# Calculating the landcover repartition in County
landcover<-data.frame(x2011 = values(nlcd_data)) #Number of rows corresponds to number of cells
landcover_freq<-table(landcover)
df_landcover <- as.data.frame(landcover_freq)
res <- df_landcover %>%
mutate(area_type_sqm = Freq*900,area_type_km=area_type_sqm/1e6,area_sqkm = sum(area_type_km))%>%
group_by(landcover)%>%
mutate(pc_land =round(100*area_type_km/area_sqkm,1))
head(arrange(res,desc(pc_land)))
# A tibble: 6 x 6
# Groups: landcover [6]
landcover Freq area_type_sqm area_type_km area_sqkm pc_land
<fct> <int> <dbl> <dbl> <dbl> <dbl>
1 52 33455938 30110344200 30110. 51107. 58.9
2 42 16073820 14466438000 14466. 51107. 28.3
3 71 5999412 5399470800 5399. 51107. 10.6
4 31 488652 439786800 440. 51107. 0.9
5 21 362722 326449800 326. 51107. 0.6
6 22 95545 85990500 86.0 51107. 0.2
## Total area calculated from raster is 51,107 square km
apache_area <- as.data.frame(apache) %>%
mutate(AREA=(as.numeric(ALAND)+as.numeric(AWATER))/1e6) %>%
select(AREA)
apache_area$AREA
29055.47 #Official area of apache county (in square km)
- 对 shapefile 和栅格进行目视检查:
差异似乎不足以证明差异是合理的
apache <- spTransform(apache,proj4string(nlcd_data))
plot(nlcd_data)
plot(apache,add=TRUE)
解决方法
原因是你得到了墨卡托投影中返回的数据。
crs(nlcd_data)
CRS arguments:
+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext
+no_defs
这个坐标参考系统保持形状,这就是它被用于网络映射的原因。它不应该用于大多数其他目的,因为它也因区域失真而臭名昭著。
实事求是的信息是永远不要只相信投影栅格的名义分辨率并假设它是正确的和/或恒定的。 可靠计算面积的方法是使用经度/纬度坐标,因为根据定义,这些坐标不会失真。
报告的空间分辨率为
res(nlcd_data)
[1] 30 30
因此,您预计单元格的面积为 30 x 30 = 900 m2 也就不足为奇了。但是阿帕奇县的像元大小实际上在 573 到 625 m2 之间。如下图。
先获取数据
library(FedData)
counties <- raster::getData("GADM",country="USA",level=2)
apache <- subset(counties,counties$NAME_2=="Apache")
nlcd <- get_nlcd(apache,year = 2011,label = "nlcd_apache",force.redo = TRUE)
# move to terra
library(terra)
r <- rast(nlcd)
ap <- vect(apache)
# county boundaries to Mercator
apm <- project(ap,crs(r))
为了计算网格单元的面积,我将它们表示为多边形。我首先聚合以获得更大的单元格以避免获得太多的小多边形(这将花费很长时间,并且可能会导致 R 崩溃)。然后我将墨卡托多边形转换为经度/纬度,计算它们的真实面积,然后转换回来(只是为了保持一致的显示目的)。
f <- 300
a <- aggregate(rast(r),f)
p <- as.polygons(a)
# compute area
g <- project(p,"+proj=longlat")
g$area <- round(expanse(g) / (f * f),5)
# project back and plot
merc <- project(g,crs(r))
plot(merc,"area",border=NA)
lines(apm)
该地图显示了原始“900 m2”像元大小的近似变化(介于 573 和 625 之间)。当我使用原始数据时,情况并非如此,如下图所示。
library(terra)
# "filename" is the local file that has the nlcd data
usnlcd <- rast(filename)
crs(usnlcd,proj4=TRUE)
#[1] "+proj=aea +lat_0=23 +lon_0=-96 +lat_1=29.5 +lat_2=45.5 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
res(x)
#[1] 30 30
请注意,+proj=aea
代表 Albers 等面积投影!
ap <- vect(apache)
apm <- project(ap,crs(usnlcd))
x <- crop(usnlcd,apm)
par(mfrow=c(1,2))
plot(x)
lines(apm)
# gg2 computed as above
plot(gg2,border=NA)
如您所见,单元面积确实是 900 m2,只有很小的失真,可以忽略不计。
您可以将墨卡托数据转换回原始的 +proj=aea
,但这样会使数据质量降低两次。所以这是一个非常糟糕的主意。您也可以考虑每个单元格的真实单元格面积,但这是相当复杂的。
最后得到每个土地覆盖类别所覆盖的面积
m <- mask(x,apm)
f <- freq(m)
f <- data.frame(f)
f$area <- round(f$count * 900 / 1e6,1)
# the next step is a bit tricky and will be done by `freq` in future versions
levs <- levels(m)[[1]]
f$class <- levs[f$value + 1]
瞧:
f[,c("class","area")]
# class area
#1 Open Water 21.5
#2 Developed,Open Space 175.3
#3 Developed,Low Intensity 46.4
#4 Developed,Medium Intensity 9.7
#5 Developed,High Intensity 1.0
#6 Barren Land 232.5
#7 Deciduous Forest 44.6
#8 Evergreen Forest 7604.6
#9 Mixed Forest 2.0
#10 Shrub/Scrub 18262.8
#11 Herbaceuous 2514.4
#12 Hay/Pasture 1.9
#13 Cultivated Crops 0.0
#14 Woody Wetlands 38.8
#15 Emergent Herbaceuous Wetlands 53.6
总数符合预期
sum(f$area)
#[1] 29009.1
附注。这个问题现在已经在源头解决了 --- 但我希望这个答案对其他使用墨卡托 CRS 空间数据的人仍然有用。
,你们大家,非常感谢您发现这个问题,也感谢 Robert 在 Github 上提醒我!我没有意识到 WCS 正在提供转换后的数据。幸运的是,MRLC 也以原始格式提供了数据,因此我已经继续向 FedData 推送更新,以提供这些数据(在 CONUS Albers 中)。
FWIW(没有图片,因为这是我的第一篇 SO 帖子):
#devtools::install_github("ropensci/FedData")
library(FedData)
library(magrittr)
library(dplyr)
library(tigris)
library(raster)
#Get Apache polygon
apache <-
tigris::counties(state = "AZ") %>%
dplyr::filter(NAME == "Apache")
# Get NCLD data
nlcd_data <-
FedData::get_nlcd(
template = apache,label = "Apache",force.redo = TRUE
)
# Transform Apache polygon to NLCD CRS
apache %<>%
sf::st_transform(raster::crs(nlcd_data))
# Plot NLCD raster and transformed Apache polygon
raster::plot(nlcd_data)
apache %>%
sf::st_geometry() %>%
plot(border = "white",lwd = 2,add = TRUE)
https://i.imgur.com/5b09t7P.png
# Mask NLCD data by Apache County (and plot result)
nlcd_data %<>%
raster::mask(apache)
plot(nlcd_data)
https://i.imgur.com/UinQ0v3.png
# Table of areas in km^2
nlcd_data_table <-
nlcd_data %>%
values() %>%
tibble::tibble(landcover = .) %>%
na.omit() %>%
dplyr::group_by(landcover) %>%
count(name = "freq") %>%
dplyr::mutate(area = (freq * 900) %>%
units::set_units("m^2") %>%
units::set_units("km^2"))
nlcd_data_table
#> # A tibble: 15 x 3
#> # Groups: landcover [15]
#> landcover freq area
#> <int> <int> [km^2]
#> 1 11 24134 21.7206
#> 2 21 194720 175.2480
#> 3 22 51321 46.1889
#> 4 23 10485 9.4365
#> 5 24 1069 0.9621
#> 6 31 259987 233.9883
#> 7 41 48906 44.0154
#> 8 42 8456196 7610.5764
#> 9 43 2179 1.9611
#> 10 52 19941185 17947.0665
#> 11 71 3190650 2871.5850
#> 12 81 2169 1.9521
#> 13 82 4 0.0036
#> 14 90 42628 38.3652
#> 15 95 59014 53.1126
# Total area NLCD raster:
sum(nlcd_data_table$area)
#> 29056.18 [km^2]
# Area of Apache county
apache %>%
sf::st_area() %>%
units::set_units("km^2")
#> 29056.19 [km^2]
由 reprex package (v2.0.0) 于 2021 年 5 月 26 日创建