从具有不同 x,y 分辨率的内核平滑创建光栅文件

Question

我有一组 GPS 数据，我正尝试使用 'kernsmooth' 包中的 bkde2d 对其进行内核平滑处理。我已经使用 'ks' 包中的 Hpi 带宽估算器来确定我的带宽，但是当我运行内核平滑并将结果列表转换为光栅时，结果产品似乎具有不同的 x、y 分辨率和因此不可能导出为 ascii。导出为 GRD 文件时也无法将此栅格读入 GIS 工具，因为它似乎已损坏，可能是由于分辨率不同。

这是我运行的一些示例代码。我的数据以 UTM30、WGS84 投影：

bnd=Hpi(x=cbind(GPS$lon, GPS$lat))
coord <- cbind(GPS$lon, GPS$lat)

est <- bkde2D(coord, bandwidth=bnd, gridsize = c(4000L, 4000L))

est.raster = raster(list(x=est$x1,y=est$x2,z=est$fhat))
projection(est.raster) <- CRS("+proj=utm +ellps=WGS84 +datum=WGS84 +zone=30 +north +units=km")`
xmin(est.raster) <- min(GPS$lon)
xmax(est.raster) <- max(GPS$lon)
ymin(est.raster) <- min(GPS$lat)
ymax(est.raster) <- max(GPS$lat)

writeRaster(est.raster, "kerntest", format='ascii')

生成的栅格图层如下所示：

class       : RasterLayer 
dimensions  : 4000, 4000, 1.6e+07  (nrow, ncol, ncell)
resolution  : 0.03242282, 0.03011303  (x, y)
extent      : 415.2883, 544.9796, 6371.946, 6492.398  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +ellps=WGS84 +datum=WGS84 +zone=30 +units=km +towgs84=0,0,0 
data source : in memory
names       : layer 
values      : 0, 0.005935748  (min, max)

然而，当我尝试导出它时，我收到错误消息：

Error in .startAsciiWriting(x, filename, ...) : 
  x has unequal horizontal and vertical resolutions. Such data cannot be stored in arc-ascii format

我的问题是为什么我的分辨率不同，我应该如何解决这个问题？

Answer 1

您可以将栅格重新采样为具有相同 x 和 y 分辨率的新栅格。您可能会以这种方式丢失一些信息。或者，您可以确保 bkde2D 中的网格大小将平分 x 和 y 范围。

est.raster <- raster::resample(est.raster, raster(ext=extent(c(415.2883, 544.9796, 6371.946, 6492.398)),resolution=0.03,crs=projection(est.raster))

Answer 2

您的示例存在几个问题（包括我们无法重现！）。它应该是这样的：

library(ks)
library(KernSmooth)
library(raster)

set.seed(0)
GPS <- data.frame(lon=runif(100), lat=runif(100)*2)
#bnd <- Hpi(GPS)
est <- bkde2D(GPS, bandwidth=0.1, gridsize = c(400L, 400L))

names(est) <- c('x', 'y', 'z')
est.raster <- raster(est)
# do not change the extent! 
projection(est.raster) <- "+proj=utm +zone=30 +north +units=km +datum=WGS84"

writeRaster(est.raster, "kerntest", format='ascii')

这个对象没有问题。但是您选择使用的文件格式无法保存这些数据。使用其他格式！例如：

writeRaster(est.raster, "kerntest.tif")

您也可以尝试使用 range.x

强制 bkde2D 生成具有方形像元的栅格

est <- bkde2D(GPS, bandwidth=0.1, gridsize = c(400L, 800L), range.x=list(c(1/800,1-1/800), c(1/800,2-1/800)))

names(est) <- c('x', 'y', 'z')
est.raster <- raster(est)
projection(est.raster) <- "+proj=utm +zone=30 +north +units=km +datum=WGS84"

writeRaster(est.raster, "kerntest", format='ascii')

从具有不同 x,y 分辨率的内核平滑创建光栅文件

Creating a raster file from a kernel smooth with differing x,y resolutions

gis

r

raster

kernel-density

r-raster