r语言 - 从 WGS84 到莫尔威德的栅格转换会改变值?

Question

我有一个关于将具有分类值的栅格从 WGS84 转换为 Mollweide 投影的问题。看起来转换会导致数据集值的更改。非常不幸的是，我很难为您提供一个可重现的示例，因此我将为您提供有关我的方法的一些详细信息。您可能有一些关于我的问题可能来自何处的提示，因为这可能是一个常见问题。欧盟网站 https://ghsl.jrc.ec.europa.eu/download.php?ds=bu 使我可以访问以下栅格：

• SMOD图层为我提供了人类住区信息(我将其转换为1(城市掩模，其中"1"表示城市地区，"NA"表示非城市地区2(农村掩码，"1"表示农村地区，"NA"表示非农村地区(。SMOD仅在Mollweide投影中可用。

• POP 层为我提供了人口密度(每个网格单元的人数(。POP 在 Mollweide 和 WGS84 中均可用。 我尝试了两种方法来估计农村和城市人口数量。

我的挑战是，对于每种方法，我都会得到不同的数字。我想知道为什么会这样：

方法1(将 SMOD 摩尔威德投影更改为 WGS84

# layers as provided by website
SMOD_MollweideProj <-  raster ("./GHS_SMOD_POP2015_GLOBE_R2019A_54009_1K_V1_0.tif") 
POP_WGSproj <- raster("./GHS_POP_E2015_GLOBE_R2019A_4326_30ss_V1_0.tif") 
SMOD_WGSproj <- projectRaster(from=SMOD_MollweideProj, to= POP_WGSproj, method='ngb' , over=T )
#create rural and urban masks - Classes 30-23-22-21 if aggregated form the "urban domain", 13-12-11-10 form the "rural  domain".
SMOD_rur_mask_1K <- SMOD_WGSproj
values(SMOD_rur_mask_1K)[values(SMOD_rur_mask_1K) >14] = NA
values(SMOD_rur_mask_1K)[values(SMOD_rur_mask_1K) <=13] =  1
SMOD_urb_mask_1K <- SMOD_WGSproj
SMOD_urb_mask_1K[SMOD_urb_mask_1K<20 ] <- NA
SMOD_urb_mask_1K[SMOD_urb_mask_1K>=21 ] <- 1
#Generate rural and urban population layers, based on total population per grid cell and rural and urban masks
POP_rur_1K_WGSproj <-  POP_WGSproj * SMOD_rur_mask_1K
POP_urb_1K_WGSproj <-  POP_WGSproj * SMOD_urb_mask_1K
#urban and rural population estimates
cellStats(POP_rur_1K_WGSproj, sum, na.rm=T) 
#2024108119 which is different to the value I get with the second approach
cellStats(POP_urb_1K_WGSproj, sum, na.rm=T)
#5321638069 which is different to the value I get with the second approach
> SMOD_MollweideProj
class      : RasterLayer 
dimensions : 18000, 36082, 649476000  (nrow, ncol, ncell)
resolution : 1000, 1000  (x, y)
extent     : -18041000, 18041000, -9e+06, 9e+06  (xmin, xmax, ymin, ymax)
crs        : +proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs 
values     : 10, 30  (min, max)
> SMOD_WGSproj
class      : RasterLayer 
dimensions : 21600, 43200, 933120000  (nrow, ncol, ncell)
resolution : 0.008333333, 0.008333333  (x, y)
extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
values     : 10, 30  (min, max)
> POP_WGSproj
class      : RasterLayer 
dimensions : 21600, 43200, 933120000  (nrow, ncol, ncell)
resolution : 0.008333333, 0.008333333  (x, y)
extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
values     : 0, 459434.6  (min, max)

方法2(使用 SMOD 和 POP Mollweide 投影

# layers as provided by website
SMOD_MollweideProj <-  raster ("./GHS_SMOD_POP2015_GLOBE_R2019A_54009_1K_V1_0.tif") 
POP_MollweideProj <- raster("./GHS_POP_E2015_GLOBE_R2019A_54009_1K_V1_0.tif") 
#create rural and urban masks - Classes 30-23-22-21 if aggregated form the "urban domain", 13-12-11-10 form the "rural  domain".
SMOD_rur_mask_1K <- SMOD_MollweideProj
values(SMOD_rur_mask_1K)[values(SMOD_rur_mask_1K) >14] = NA
values(SMOD_rur_mask_1K)[values(SMOD_rur_mask_1K) <=13] =  1
SMOD_urb_mask_1K <- SMOD_MollweideProj
SMOD_urb_mask_1K[SMOD_urb_mask_1K<20 ] <- NA
SMOD_urb_mask_1K[SMOD_urb_mask_1K>=21 ] <- 1
#Generate rural and urban population layers, based on total population per grid cell and rural and urban masks
POP_rur_1K_MollweideProj <-  POP_MollweideProj * SMOD_rur_mask_1K
POP_urb_1K_MollweideProj <-  POP_MollweideProj * SMOD_urb_mask_1K
#urban and rural population estimates
cellStats(POP_rur_1K_MollweideProj, sum, na.rm=T)
# 1726372189 which is different to the value I get with the first approach
cellStats(POP_urb_1K_MollweideProj, sum, na.rm=T)
# 5622956252 which is different to the value I get with the first approach

非常感谢您的建议

Answer 1

根据上面的评论，这是我为您的方法 2 建议的代码，即使用 SMOD 和 POP 的 Mollweide 投影。我下载了单个单元格而不是全局层的数据，只是为了减少执行时间。

请注意，我用raster::mask()将 POP 屏蔽到城市和农村。使用此方法，无需在蒙版中设置1值，您可以在设置为NA要遮罩的单元格后保留原始值。请参阅?raster::mask。

library(raster)
smod <- raster("data/GHS_SMOD_POP2015_GLOBE_R2019A_54009_1K_V1_0_11_4.tif")
pop <- raster("data/GHS_POP_E2015_GLOBE_R2019A_54009_1K_V1_0_11_4.tif")
# produce rural mask
maskRural <- smod
maskRural[maskRural > 14] <- NA
# produce urban mask
maskUrban <- smod
maskUrban[maskUrban < 20] <- NA
# use raster::mask to produce rural and urban population layers
popRural <- raster::mask(pop, maskRural)
popUrban <- raster::mask(pop, maskUrban)
# check that the total population of rural + urban is equal to the original
sum(pop[], na.rm = TRUE)
sum(popUrban[], na.rm = TRUE) + sum(popRural[], na.rm = TRUE)

最后几行只是检查农村和城市人口的总和是否等于原始总人口。因此，对于我使用的单个单元格，结果是：

> sum(pop[], na.rm = TRUE)
[1] 67487835
> sum(popUrban[], na.rm = TRUE) + sum(popRural[], na.rm = TRUE)
[1] 67487835

要投影到 WGS84，您可以执行以下操作：

popUrbanWgs84 <- projectRaster(popUrban, crs = crs("+init=epsg:4326"), method = "bilinear")

您也可以在此处指定一个res参数，否则它将为您选择一个，但问题是使用什么分辨率？1公里的原始分辨率约为0.008度，但就经度而言，这在全球范围内有所不同。

我的建议是，如果可能的话，坚持使用Mollweide。