根据 NLCD（土地覆盖数据）栅格计算的县域太大答案

【问题标题】：County area calculated from NLCD (Landcover data) rasters is too large根据 NLCD（土地覆盖数据）栅格计算的县域太大
【发布时间】：2021-08-14 15:26:58
【问题描述】：

我正在尝试计算美国每个县的土地覆盖重新分配。我已经使用 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)

【问题讨论】：

标签： r gis raster area

【解决方案1】：

原因是你得到了墨卡托投影返回的数据。

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 m² 也就不足为奇了。但阿帕奇县的单元大小实际上在 573 到 625 m² 之间。如下图所示。

先获取数据

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 m²”单元（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 代表阿尔伯斯 等面积投影！

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, "area", border=NA)

如您所见，单元格面积确实是900 m²，失真很小，小到可以忽略不计。

您可以将墨卡托数据转换回原始的+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

PS。这个问题现在已经从源头上解决了 --- 但我希望这个答案对使用墨卡托 CRS 空间数据的其他人仍然有用。

【讨论】：

非常感谢您的详细解答。我不熟悉 terra 库，所以这段代码目前对我来说是神秘的，但我会完成它。您能否确认我理解这一点：实现我想要做的更好的方法是使用原始 nlcd 数据，根据原始 nlcd 数据的 crs 转换我的 shapefile，裁剪县并进行计算从那里？
是的。在crop 之后，您可以使用mask(x, apm) 将Apache 之外的单元格设置为NA。

【解决方案2】：

你们，非常感谢你发现这个，感谢 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,
    year = 2011,
    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-05-26 创建}

【讨论】：