重新投影光栅数据时极值较小

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”功能的数据。这一定是一个常见的问题，但我在网上找不到任何帮助。

非常感谢。

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