I have a set of GPS data which I am attempting to kernel smooth using the bkde2d in the 'kernsmooth' package. I have used the Hpi bandwidth estimator in the 'ks' package to determine my bandwidth however when I run the kernel smooth and convert the resulting list into a raster the resulting product appears to have differing x, y resolutions and is therefore impossible to export as an ascii. It is also impossible to read this raster into a GIS tool when exporting as a GRD file as it appears to be corrupt, presumably due to having differing resolutions.
Here is some sample code from my run. My data is projected in 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')
The resulting raster layer looks like this:
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)
However when I attempt to export it I get the error message:
Error in .startAsciiWriting(x, filename, ...) :
x has unequal horizontal and vertical resolutions. Such data cannot be stored in arc-ascii format
My question is why are my resolutions different and how should I resolve this issue?
You can resample the raster to a new raster with equal x and y resolutions. You may lose some information this way. Alternatively you could make sure your gridsize in bkde2D would divide the x and y extents equally.
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))