从 WGS84 到 Mollweide 的光栅转换会改变值吗?
raster conversion from WGS84 to Mollweide alters values?
我有一个关于将具有分类值的栅格从 WGS84 转换为 Mollweide 投影的问题。看起来转换会导致数据集值的改变。非常遗憾,我很难为您提供可重现的示例,因此我将为您提供有关我的方法的一些细节。您可能对我的问题可能来自何处有一些提示,因为这可能是一个常见问题。欧盟网站 https://ghsl.jrc.ec.europa.eu/download.php?ds=bu 允许我访问以下栅格:
SMOD 图层为我提供了人类住区信息(我将 SMOD 数据集转换为 1)城市掩码,其中“1”表示城市地区,“NA”表示非城市地区 2)农村掩码农村地区为“1”,非农村地区为“NA”)。 SMOD 仅适用于 Mollweide 投影。
POP 层为我提供了人口密度(每个网格单元的人数)。 Mollweide 和 WGS84 都提供 POP。
我尝试了两种方法来估计农村和城市人口数量。
我面临的挑战是,我对这些方法中的每一种都得到了不同的数字。我想知道为什么会这样:
方法 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。
我有一个关于将具有分类值的栅格从 WGS84 转换为 Mollweide 投影的问题。看起来转换会导致数据集值的改变。非常遗憾,我很难为您提供可重现的示例,因此我将为您提供有关我的方法的一些细节。您可能对我的问题可能来自何处有一些提示,因为这可能是一个常见问题。欧盟网站 https://ghsl.jrc.ec.europa.eu/download.php?ds=bu 允许我访问以下栅格:
SMOD 图层为我提供了人类住区信息(我将 SMOD 数据集转换为 1)城市掩码,其中“1”表示城市地区,“NA”表示非城市地区 2)农村掩码农村地区为“1”,非农村地区为“NA”)。 SMOD 仅适用于 Mollweide 投影。
POP 层为我提供了人口密度(每个网格单元的人数)。 Mollweide 和 WGS84 都提供 POP。 我尝试了两种方法来估计农村和城市人口数量。
我面临的挑战是,我对这些方法中的每一种都得到了不同的数字。我想知道为什么会这样:
方法 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。