从 NLCD（土地覆盖数据）栅格计算的国家面积太大
County area calculated from NLCD (Landcover data) rasters is too large
我正在尝试计算美国每个县的土地覆被重新划分。
我已经使用 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
- 然后我提取了栅格的值并将它们放入频率 table 中。从那里我计算结果面积。由于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
此坐标参考系统保留形状，这就是它用于网络制图的原因。它不应该用于大多数其他目的，因为它也因 区域变形 .
而臭名昭著。
要点是永远不要只相信投影光栅的标称分辨率并假设它是正确的 and/or 常数。 可靠 计算面积的方法是使用 longitude/latitude 坐标，因为根据定义，这些坐标没有变形。
报告的空间分辨率为
res(nlcd_data)
[1] 30 30
因此，您预期单元格的面积为 30 x 30 = 900 m2 也就不足为奇了。但 Apache 县的像元大小实际上在 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 崩溃）。然后我将墨卡托多边形转换为 longitude/latitude，计算它们的真实面积，然后转换回来（只是为了一致的显示目的）。
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, "area", border=NA)
可以看到，cell面积确实有900m2，只有很小的畸变，小到可以忽略
您可以将墨卡托数据转换回原始数据 +proj=aea，但这样会使数据质量降低两次。所以这是一个非常糟糕的主意。您可以 还可以计算每个单元格的真实单元格面积，但这相当复杂。
最后，得到每个土地覆盖的面积class
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 空间数据的人仍然有用。
你们好，非常感谢你看到这个，也感谢 Robert 在 Github 上提醒我！我没有意识到 WCS 正在提供转换后的数据。幸运的是，MRLC 也以其原始格式提供了数据，因此我继续推进并推送了对提供这些数据的 FedData 的更新（在 CONUS Albers 中）。
FWIW（没有图片，因为这是我的第一个 SO post）：
#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 创建
我正在尝试计算美国每个县的土地覆被重新划分。 我已经使用 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
- 然后我提取了栅格的值并将它们放入频率 table 中。从那里我计算结果面积。由于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
此坐标参考系统保留形状，这就是它用于网络制图的原因。它不应该用于大多数其他目的，因为它也因 区域变形 .
而臭名昭著。要点是永远不要只相信投影光栅的标称分辨率并假设它是正确的 and/or 常数。 可靠 计算面积的方法是使用 longitude/latitude 坐标，因为根据定义，这些坐标没有变形。
报告的空间分辨率为
res(nlcd_data)
[1] 30 30
因此，您预期单元格的面积为 30 x 30 = 900 m2 也就不足为奇了。但 Apache 县的像元大小实际上在 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 崩溃）。然后我将墨卡托多边形转换为 longitude/latitude，计算它们的真实面积，然后转换回来（只是为了一致的显示目的）。
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, "area", border=NA)
可以看到，cell面积确实有900m2，只有很小的畸变，小到可以忽略
您可以将墨卡托数据转换回原始数据 +proj=aea，但这样会使数据质量降低两次。所以这是一个非常糟糕的主意。您可以 还可以计算每个单元格的真实单元格面积，但这相当复杂。
最后，得到每个土地覆盖的面积class
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 空间数据的人仍然有用。
你们好，非常感谢你看到这个，也感谢 Robert 在 Github 上提醒我！我没有意识到 WCS 正在提供转换后的数据。幸运的是，MRLC 也以其原始格式提供了数据，因此我继续推进并推送了对提供这些数据的 FedData 的更新（在 CONUS Albers 中）。
FWIW（没有图片，因为这是我的第一个 SO post）：
#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 创建