Search code examples
rspatialraster

R Assigning x or y coordinate to cells of a raster to perform calculations


Is there any solution to assign to assign X or Y coordinate to all of the cells of a raster image using R?

For example, suppose I have a raster with 3x3 cells. The lower left coordinate is X=7,Y=15 (in meters), and resolution is +10 m for both, X and Y direction (X increases rightward, and Y upward). Then, I would like to generate raster surfaces where each cell has the X and Y value, like these ones:

X raster
7  17  27
7  17  27
7  17  27

Y surface
35  35  35
25  25  25
15  15  15

UPDATE: this is the actual raster object.

Is there any way to do this?

I tried the package 'raster', but couldn't find a solution.

Any help would be appreciated.


Solution

  • This is the solution that I found, inspired in the comment submitted by @Carl.

    Suppose my raster is called d, with these characteristics:

    class       : RasterLayer
    dimensions  : 59, 67, 3953  (nrow, ncol, ncell)
    resolution  : 90, 90  (x, y)
    extent      : 482855.6, 488885.6, 4763517, 4768827  (xmin, xmax, ymin, ymax)
    coord. ref. : +proj=tmerc +lat_0=0 +lon_0=-3 +k=1 +x_0=500000 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
    data source : in memory
    names       : subset
    values      : 328, 1120  (min, max)
    

    I first generated a matrix with the x coordinates values:

    xm<-matrix(xFromCell(d,c(1:3953)),nrow=59,byrow=TRUE)
    

    Then created a raster with the matrix:

    x<-raster(xm,xmn=482855.6, xmx=488885.6,ymn=4763517,ymx=4768827)
    

    And finally, assigned its projection:

    projection(x)<-"+proj=tmerc +lat_0=0 +lon_0=-3 +k=1 +x_0=500000 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
    

    I also displayed the raster, to see if results where OK...

    plot(x)
    

    raster image where each cell has the x-coordinate value

    ...and read the raster description

    x
    class       : RasterLayer
    dimensions  : 59, 67, 3953  (nrow, ncol, ncell)
    resolution  : 90, 90  (x, y)
    extent      : 482855.6, 488885.6, 4763517, 4768827  (xmin, xmax, ymin, ymax)
    coord. ref. : +proj=tmerc +lat_0=0 +lon_0=-3 +k=1 +x_0=500000 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
    data source : in memory
    names       : layer
    values      : 482900.6, 488840.6  (min, max)
    

    I repeated these steps, but using yFromCell in the first step.