问题描述
我正在尝试以0.0417度的分辨率重新投影从经度/纬度到Behrmanns等面积（epsg：6933）的人工发光的整体栅格。由于在重新投影过程中对像素进行插值时，市区周围的数据激增，整个图层上的数据丢失率约为15％。
我尝试将栅格转换为空间点数据框，重新投影空间点数据框，然后使用使用“投影仪”功能作为模板栅格创建的栅格进行栅格化（我认为模板栅格的尺寸，范围和分辨率可能是问题吗？）但是，这会产生一个带有水平线穿过图层的栅格。
以下是一些示例代码，以西班牙为例。我可以通过电子邮件将西班牙的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”函数的数据。这肯定是一个常见的问题，但我在网上找不到任何帮助。
非常感谢。
Example of horizontal lines after reprojecting stack.imgur.com/IV0fZ.png
解决方法
变换栅格时，将计算新的像元值。通常，这是通过平均来完成的，这会导致较小的极值。
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,3))
r <- round(focal(r,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,ncell)
#resolution : 482000,638000 (x,y)
#extent : -18813530,18782470,-8607770,8618230 (xmin,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,max)
或者确实使用多边形（如果您没有太多的单元格）
p <- as(r,"SpatialPolygonsDataFrame")
y <- spTransform(p,behrmann)
y
#class : SpatialPolygonsDataFrame
#features : 2160
#extent : -17367530,17367530,-7089914,7089914 (xmin,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