Search code examples
rrasterr-sf

Convert between individual STARS pixels and SF lat/long coordinates in R


I have a STARS raster object raster in R (e.g of size 100x100) and can get the coordinates using st_coordinates(raster)

But I'm having problems converting between the two. For example, I'd like to know the lat/long coordinates for only a particular raster pixel, e.g. (50, 50).

I was assuming that st_coordinates would give me a 1D array, so I can simply convert from the 2D raster matrix to the 1D array (e.g. by converting the 2D index (50, 50) into a 1D index using something like #columns*i+j, which in my example is 100*50+50).

This isn't working though, which makes me think I'm misunderstanding how the STARS raster object and coordinates indices map onto each other. Does anyone know how to bridge between the two, on a pixel-by-coordinate basis?

===

Updated with an example of what I'm trying to do:

    r <- read_stars(system.file("tif/L7_ETMs.tif", package = "stars"))

    ## Let's say I want to grab the pixel 50x50 from the first band. 
    ## I can do that going into the raster matrix like this:

    pixel <- r[[1]][50,50,1]

   ## now I can get the coordinates for the raster
   ## using sf like this
   ## this gives me a coordinate df with columns x, y and band.
 
    coordinates <- st_coordinates(r)

   ## what I want to do now is get the coordinate for the
   ## pixel I listed earlier, (50, 50, 1)  -- in my case, I only
   ## have one band so I'm not going to worry about band index logic

   ## to do this, I assume that (50, 50) is the same as
   ## ncol(r)*(50-1)+50 or 17298. I say 50-1 since R is using 1 indexing.
   ## I use this to get a coordinate from coordinates like this:

   coord <- coordinates[ncol(r)*(50-1)+50,]

   ## This returns me the following:
   

>    17298 294376.5 9119350    1

   ## If I wanted to do this many times, I could make a list 
   ## of coordinates for various pixels, then put them into 
   ## a new sf object using st_as_sf(...)

When I tried doing the above in a loop and plotting the results, there was a substantial mismatch... the raster pixels did not map to the right coordinates after plotting them in a new sf object. This has me thinking my conversion from raster 2D array to coordinate 1D list is not right. In fact, I'm realizing I have no idea at all what logic sf uses to convert the raster to a 1D list, which might explain the problem... Do you have any ideas about how these map to each other and how to index the coordinates array for a given raster pixel? Please let me know if I still need to clarify further. Thanks!


Solution

  • I think the key here is to subset a stars object then feed the subsetted object into st_coordinates.

    r <- read_stars(system.file("tif/L7_ETMs.tif", package = "stars"))
    
    # extract band 1 raster value at position 50, 50
    # two methods with same result
    
    # convert to matrix then subset
    r[[1]][50,50,1] 
    #[1] 56
    
    
    # subset stars raster then extract value
    r[,50,50,1][[1]] 
    #, , 1
    #
    #     [,1]
    #[1,]   56
    

    If you use the subset stars raster then extract value workflow, you can use the subsetted raster in the st_coordinates function.

    st_coordinates(r[,50,50,1])
    #        x       y band
    #1 290187 9119350    1