从 WGS84 到 Mollweide 的光栅转换会改变值吗?

raster conversion from WGS84 to Mollweide alters values?

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

我面临的挑战是,我对这些方法中的每一种都得到了不同的数字。我想知道为什么会这样:

方法 1) 将 SMOD Mollweide 投影更改为 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

非常感谢您的建议

根据上面的评论,这是我为您的方法 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。