How to make RASTER from irregular point data without interpolation

Question

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?

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, Nov 30 '13 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, Nov 13 '14 at 10:48

Robert Hijmans · Accepted Answer · 2023-03-17T18:36:59.907

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)

edited Mar 17 '23 at 18:36

answered Dec 01 '13 at 05:20

Robert Hijmans

10,683
25
35

excellent solution to the problem@Robert – ToNoY Dec 01 '13 at 16:16
Sorry can you clarify why you add e <- e + 1000 ? – mmann1123 Sep 10 '15 at 03: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. – Robert Hijmans Sep 10 '15 at 03:12
Is there a pythonic solution for this? – raaj Nov 12 '19 at 11:37
why in raster(e, ncol=10, nrow=2) ncol is 10 and nrow is 2? – yuliaUU May 06 '22 at 05:59
@yuliaUU: It is just an example – Robert Hijmans May 12 '22 at 00:47
@RobertHijmans how can we do this in terra? – Herman Toothrot Mar 17 '23 at 16:37
1

Now updated for terra – Robert Hijmans Mar 18 '23 at 05:53
@RobertHijmans thanks a lot. Does rasterize automatically ignore NAs? I don't see na.rm in the function options. – Herman Toothrot May 15 '23 at 16:48
No, but you can use it because there is argument "... additional arguments passed to fun" – Robert Hijmans May 15 '23 at 19:37

Tomas · Answer 2 · 2019-08-21T14:16:03.477

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.

edited Aug 21 '19 at 14:16

answered Nov 13 '14 at 10:56

Tomas

1,252
3
17
40

1

The first line of code has a closing bracket missing? – Antti Aug 21 '19 at 08:57
@Antti thanks, corrected! – Tomas Aug 21 '19 at 14:16

How to make RASTER from irregular point data without interpolation

2 Answers2

you need to provide a function 'fun' for when there are multiple points per cell

or r <- raster(xmn=, xmx=, ...

you need to provide a function 'fun' for when there are multiple points per cell

Linked