2

I was doing a query to extract pixel value from a raster and came up to some strange results. If I extracted the value using R and the Raster package, the value was different from the same point extracted using PostGIS.

Doing some further investigation, I found out that the differences only occurred when I used 'round' coordinates, that is -55.0, -22.0. These coordinates are in between two columns in my raster. So what I figured was that R was choosing the pixel from the right side of the edge and PostGIS, from the left.

Does each software round the coordinates in a different way?

7

Yes, it will be different in different software. For R the best way to understand is to experiment and dig under the hood:

r <- setExtent(raster(matrix(1:12, 3, 4)), extent(0, 4, 0, 3))
plot(r)
extract(r, cbind(4, 3))
#[1] 10 
extract(r, cbind(3, 3))
#[1] 10 
extract(r, cbind(2.99, 3))
#[1] 7

It also changes over time while software is developed, I'm pretty sure that the edge value at 4,3 used to return NA (for example).

The crux difference will usually be about whether the cell is treated as a discrete unit with exact extents (like here, and like in GDAL's default), or as a centre point with no exacting boundary, where the raster is treated like a field varying continuously.

You can see a raster treated this way in R though, try extract(..., method = "bilinear") in the examples above. It depends what the operation is, and you'll find this also varies in other software. It's actually not clearly defined if you look at the wider variety of conventions in how data is dealt with, though this "cell as rectangle" convention is pretty solid in GDAL and R, it falls down when encountering formats that store each axis-pixel's coordinate (which is redundant) or when cell-edge or cell-centre alignment is ambiguous, or just incorrect. It happens often in the more easily-created formats like NetCDF.

To see how 'raster' does it, run down the rabbit hole (type each line, it will lead you to type the next line in turn):

extract

showMethods("extract")

findMethods("extract")[["Raster#matrix"]]

raster:::.xyValues

cellFromXY

raster:::.cellValues

raster:::.doCellFromXY ... whoops now it's C++

https://github.com/cran/raster/blob/master/src/xyCell.cpp#L8

You can do it yourself using arithmetic with coordinates compared to the x and y axis positions - it's worth doing it from scratch for yourself - the key is what cell would the difference in my position from the start edge put me in, and what row you are in (or col depending on how you orient it). You'll find similar code at the bottom of every software implementation, with variations on how they store the "transform", basically the information that encodes the extent, number of pixels, and the size of each pixel.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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