重投影栅格数据时的极值较少

Question

我正在尝试以 0.0417 度的分辨率将人工光发射的全球栅格从经度/纬度重新投影到贝尔曼等面积 (EPSG:6933)。由于在重投影期间对像素进行插值时，城市区域周围的数据会出现峰值，因此整个图层会丢失大约 15% 的数据。

我尝试将栅格转换为空间点数据框，重新投影空间点数据框，然后使用使用“projectraster”函数创建的栅格作为模板栅格进行栅格化（我认为模板栅格的尺寸、范围和分辨率可能是问题所在？）但是，这会产生一个水平线穿过图层的栅格。

这里是一些以西班牙为例的示例代码。我可以通过电子邮件发送西班牙的 tif 文件 (246kb)：

library("sf")
library("raster")

behrmann <- CRS('+proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +datum=WGS84 +ellps=WGS84 +units=m +no_defs')

r <- raster("~/Documents/R spatial data[enter image description here][1]/Spain.tif")
cellStats(r, sum) # check summed light emissions
r_temp <- projectRaster(r, crs = behrmann) # creates template for rasterisation (data is lost due to interpolation of data spikes)
spdf <- rasterToPoints(r, spatial = TRUE)
spdf2 <- spTransform(spdf, CRS = behrmann)
r2 <- rasterize(spdf2, r_temp, field = "Spain", fun = "sum")
cellStats(r2, sum) # check no data has been lost
plot(log10(r2)) # see attached image[enter image description here][1]

如何在不丢失数据和避免水平线的情况下重新投影到相等的区域？我还尝试转换为空间多边形数据框而不是空间点，这不会产生线条，而是会丢失类似于“projectRaster”函数的数据。这一定是一个常见的（ish）问题，但我在网上找不到任何帮助。

非常感谢。

Example of horizontal lines after reprojectingstack.imgur.com/IV0fZ.png

Answer 1

转换栅格时，会计算新的像元值。这通常是通过平均来完成的，这样可以减少极值。

library(raster)
r <- raster(res=5)
set.seed(1)
values(r) <- runif(ncell(r))
r <- focal(r, w=matrix(1, 3, 3))
r <- focal(r, w=matrix(1, 3, 3))
r <- round(focal(r, w=matrix(1, 3, 3)))
r
#class      : RasterLayer 
#dimensions : 36, 72, 2592  (nrow, ncol, ncell)
#resolution : 5, 5  (x, y)
#extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
#crs        : +proj=longlat +datum=WGS84 +no_defs 
#source     : memory
#names      : layer 
#values     : 229, 473  (min, max)


behrmann <- "+proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +datum=WGS84 +ellps=WGS84 +units=m"
x <- projectRaster(r, crs=behrmann)
x
#class      : RasterLayer 
#dimensions : 27, 78, 2106  (nrow, ncol, ncell)
#resolution : 482000, 638000  (x, y)
#extent     : -18813530, 18782470, -8607770, 8618230  (xmin, xmax, ymin, ymax)
#crs        : +proj=cea +lat_ts=30 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
#source     : memory
#names      : layer 
#values     : 233.0725, 471.7214  (min, max)

但是，您可以改用最近邻方法（这类似于您通过转换点数据尝试的方法）。

z <- projectRaster(r, crs=behrmann, method="ngb")
z
#class      : RasterLayer 
# ...
#values     : 229, 473  (min, max)

或者确实使用多边形（如果您没有太多单元格）

p <- as(r, "SpatialPolygonsDataFrame")
y <- spTransform(p, behrmann)

y
#class       : SpatialPolygonsDataFrame 
#features    : 2160 
#extent      : -17367530, 17367530, -7089914, 7089914  (xmin, xmax, ymin, ymax)
#crs         : +proj=cea +lat_ts=30 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
#variables   : 1
#names       : layer 
#min values  :   229 
#max values  :   473