20

I was trying to make a raster image from an irregularly spaced point database. The data looks like-

> head(s100_ras)
         x       y         z
1 267573.9 2633781 213.29545
2 262224.4 2633781  69.78261
3 263742.7 2633781  51.21951
4 259328.4 2633781 301.98413
5 264109.8 2633781 141.72414
6 255094.8 2633781  88.90244

I want these 'z' values within a mesh which I created by

# Create a fine mesh grid
my_mesh=expand.grid(seq(min(s100_ras$Y),max(s100_ras$Y),l=100),
                    seq(min(s100_ras$X),max(s100_ras$X),l=100))

I also want the z-values to be assigned as 'NA' for those mesh points that are outside the data points. The points over the mesh looks like this: https://drive.google.com/file/d/0B6GUNg-8d30vYzlwTkhvaHBFTnc/edit?usp=sharing when I plot

plot(my_mesh)
points(s100_ras$Y, s100_ras$X, pch="*", col='blue')

The problem is that I'm not sure how to build on this, the following steps don't work because my mesh grid and data points are not of the same scale!!

library(rgdal)
library(raster)
xyz<-cbind(my_mesh, s100_ras)
r <- rasterFromXYZ(xyz)
image(r)

If I try to make a raster just by using the data points (without any mesh), R throws an error since my data is irregularly spaced!

library(sp)
s100_ras <- data.frame(expand.grid(x = s100_ras$Y, y = s100_ras$X), 
                       z = as.vector(s100_ras$mean))
coordinates(s100_ras) <- ~x+y
proj4string(s100_ras) <- CRS("+proj=utm +zone=46 +datum=WGS84")
gridded(s100_ras) = TRUE

suggested tolerance minimum: 0.916421 
Error in points2grid(points, tolerance, round) : 
  dimension 1 : coordinate intervals are not constant

Moreover, I was trying to play with 'rasterize' function (for irregular grids) of 'raster package', but couldn't get a way with it :(. I know how to interpolate and make a regular grid, but for the sake of originality, I want to AVOID interpolation. Is it possible to make a raster of irregularly spaced data points without idw or kriging methods?

3
  • The problem with irregular spaced grids is that the algorithm fails if points lie too close/far together. A (non-optimal) workaround: Why not take the minimal distance between cells and create a rectangular vector grid. Then join the average values of the points to that grid and rasterize it.
    – Curlew
    Commented Nov 30, 2013 at 23:20
  • I had the same problem - the solution was to use SpatialPixelsDataFrame with the suggested tolerance argument (0.916421 in your case).
    – Tomas
    Commented Nov 13, 2014 at 10:48
  • @Tomas how do you decide on the tolerance? Commented May 15, 2023 at 16:35

2 Answers 2

26

I assume you want your irregular point data on a regular raster. In that case, rasterize should work, and the examples in ?rasterize show how.

Your example data

s100 <- matrix(c(267573.9, 2633781, 213.29545, 262224.4, 2633781, 69.78261, 263742.7, 2633781, 51.21951, 259328.4, 2633781, 301.98413, 264109.8, 2633781, 141.72414, 255094.8, 2633781, 88.90244),  ncol=3,  byrow=TRUE)
colnames(s100) <- c('X', 'Y', 'Z')

Using the "terra" package:

library(terra)
# set up a raster geometry, here deriving an extent from your data
e <- ext(apply(s100[,1:2], 2, range))
# set up the raster, for example
r <- rast(e, ncol=10, nrow=2)

# you need to provide a function 'fun' for when there are multiple points per cell
x <- rasterize(s100[, 1:2], r, s100[,3], fun=mean)
plot(x)

And with the old raster package

library(raster)
# set up an 'empty' raster, here via an extent object derived from your data
e <- extent(s100[,1:2])
e <- e + 1000 # add this as all y's are the same

r <- raster(e, ncol=10, nrow=2)
# or r <- raster(xmn=, xmx=,  ...

# you need to provide a function 'fun' for when there are multiple points per cell
x <- rasterize(s100[, 1:2], r, s100[,3], fun=mean)
plot(x)
10
  • excellent solution to the problem@Robert
    – ToNoY
    Commented Dec 1, 2013 at 16:16
  • Sorry can you clarify why you add e <- e + 1000 ?
    – mmann1123
    Commented Sep 10, 2015 at 3:00
  • @mman1123 That is just to make things work with these odd example data. All y coordinates are the same and hence the extent is zero in the y direction, so I add 1000 --- entirely arbitrary --- to be able to make a raster from the extent. Commented Sep 10, 2015 at 3:12
  • Is there a pythonic solution for this?
    – raaj
    Commented Nov 12, 2019 at 11:37
  • 1
    Now updated for terra Commented Mar 18, 2023 at 5:53
6

This worked for me - the solution was to use SpatialPixelsDataFrame with the suggested tolerance argument (0.916421 in your case):

points <- SpatialPoints(s100_ras[,c('x','y')], s100_ras[,c('z')])
pixels <- SpatialPixelsDataFrame(points, tolerance = 0.916421, points@data)
raster <- raster(pixels[,'z'])

though, due to the high tolerance value, the raster doesn't fit very well to the original points. Could be fit much better.

2
  • 1
    The first line of code has a closing bracket missing?
    – Antti
    Commented Aug 21, 2019 at 8:57
  • @Antti thanks, corrected!
    – Tomas
    Commented Aug 21, 2019 at 14:16

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.