library(raster)
library(rnaturalearth)
library(terra)
r <- raster::getData('CMIP5', var='tmin', res=10, rcp=45, model='HE', year=70)
r <- r[[1]]
shp <- rnaturalearth::ne_countries()
newcrs <- "+proj=robin +datum=WGS84"
r <- rast(r)
shp <- vect(shp)
r_pr <- terra::project(r, newcrs)
shp_pr <- terra::project(shp, newcrs)
对于shp_pr
中的每个国家,我想标准化底层栅格0-1分。这意味着将一个单元格除以一个国家边界内所有单元格的总和,并对所有国家重复此方法。我的操作如下:
country_vec <- shp$sovereignt
temp_ls <- list()
for(c in seq_along(country_vec)){
country_ref <- country_vec[c]
if(country_ref == "Antarctica") { next }
shp_ct <- shp[shp$sovereignt == country_ref]
r_country <- terra::crop(r, shp_ct) # crops to the extent of boundary
r_country <- terra::extract(r_country, shp_ct, xy=T)
r_country$score_norm <- r_country$he45tn701/sum(na.omit(r_country$he45tn701))
r_country_norm_rast <- rasterFromXYZ(r_country[ , c("x","y","score_norm")])
temp_ls[[c]] <- r_country_norm_rast
rm(shp_ct, r_country, r_country_norm_rast)
}
m <- do.call(merge, temp_ls)
我想知道这是否是最有效/正确的方法来做到这一点,即没有任何for循环,有人有任何建议吗?
稍微更新和简化的示例数据(不需要投影数据)
library(terra)
library(geodata)
r <- geodata::cmip6_world("HadGEM3-GC31-LL", "585", "2061-2080", "tmin", 10, ".")[[1]]
v <- world(path=".")
v$ID <- 1:nrow(v)
解决方案
z <- rasterize(v, r, "ID", touches=TRUE)
zmin <- zonal(r, z, min, na.rm=TRUE, as.raster=TRUE)
zmax <- zonal(r, z, max, na.rm=TRUE, as.raster=TRUE)
x <- (r - zmin) / (zmax - zmin)
请注意,上面的代码对0到1之间的每个国家的单元格值进行了标准化。
要转换数据,使值的总和为1(按国家),您可以执行:
z <- rasterize(v, r, "ID", touches=TRUE)
zsum <- zonal(r, z, sum, na.rm=TRUE, as.raster=TRUE)
x <- r / zsum