Skip to content

Extract xy coordinates from raster of interesting neighborhood cells

4 messages · ASANTOS, Ben Tupper, Vijay Lulla

Original
ASANTOS
# 
Dear R-Sig-Geo Members,

 ??? I've like to extract xy coordinates from raster of 24 neigborhood 
cells that surround the 4 given points (pts_s) pixels and for this I used:

#Packages

require(raster)
require(sp)

## Create a raster
r <- raster(nc=30, nr=30)
r <- setValues(r, round(runif(ncell(r))* 255))
plot(r)

##Given interesting points coordinates
xd???? <- c(-24.99270,45.12069,99.40321,73.64419)
yd? <- c(-45.435267,-88.369745,-7.086949,44.174530)
pts <- data.frame(xd,yd)
pts_s<- SpatialPoints(pts)
points(pts_s, col="red", pch=16)

## Find pixels center of each point (My idea is not to use the original 
point coordinates, but coordinates of the center of interesting pixel).
N_cells <- cellFromXY(r, pts_s)


# Extract xy coordinates from raster of 24 neighborhood cells given 
pixel number N_cells
e<-adjacent(r, N_cells , directions='bishop', id=TRUE)
coordinates(land_mask)[e[,1],] ## Doesn't return 24 neighborhood by 
N_cells object
ng_coords<-coordinates(land_mask)[e[,1],]
#

#Visualization

plot(r)
points(pts_s, col="red", pch=16)
points(ng_coords, col="black", pch=16)
##

 ??? But I don't have success because the ng_coords object is not a 24 
neighborhood cells of each point that I search. There are solution for this?

Thanks in advance,

Alexandre
Ben Tupper
# 
Hi,

I'm not sure why you expect 24 points - you have 4 locations and for each you want the 4 bishop's directions - so, I think at most you should expect 16 points.  See http://www.lpc.uottawa.ca/publications/moransi/moran.htm <http://www.lpc.uottawa.ca/publications/moransi/moran.htm> for where I expect 4 points for bishop's for each input location; perhaps you have a different idea of bishop's direction?  

Your second point [45.12069, -88.369745] is on to the southern edge of the raster - so it can't find points southward - only the 2 northward.

Finally, it may not make any difference at all, but your pts_s has no coordinate reference.  Below I show where I assign the same projection as your raster.  I also don't have 'land_mask' so I reused 'r' and used xyFromCell to backwards extract the adjacent cell coordinates.

Does that do the trick for you?

Cheers,
Ben

### START
require(raster)
require(sp)

## Create a raster
r <- raster(nc=30, nr=30)
r <- setValues(r, round(runif(ncell(r))* 255))
plot(r)

##Given interesting points coordinates
xd     <- c( -24.99270,   45.12069,  99.40321,  73.64419)
yd     <- c(-45.435267, -88.369745, -7.086949, 44.174530)
pts <- data.frame(xd,yd)
pts_s<- SpatialPoints(pts)
projection(pts_s) <- projection(r)
points(pts_s, col="red", pch=16)

## Find pixels center of each point 
N_cells <- cellFromXY(r, pts_s)

e <- adjacent(r, N_cells , directions='bishop', id=TRUE, sorted = TRUE)
xy <- xyFromCell(r, e[,'to'], spatial = TRUE)


#Visualization

plot(r)
points(pts_s, col="red", pch=16)
points(xy, col="black", pch=16)
### END
Ben Tupper
Bigelow Laboratory for Ocean Sciences
60 Bigelow Drive, P.O. Box 380
East Boothbay, Maine 04544
http://www.bigelow.org

Ecological Forecasting: https://eco.bigelow.org/
Vijay Lulla
# 
Or you can try something like this?

### START
require(raster)
require(sp)

r <- raster(nc=30, nr=30)
r <- setValues(r, round(runif(ncell(r))* 255))

xd     <- c(-24.99270,45.12069,99.40321,73.64419)
yd  <- c(-45.435267,-88.369745,-7.086949,44.174530)
pts <- data.frame(xd,yd)
pts_s<- SpatialPoints(pts)

N_cells <- cellFromXY(r, pts_s)
e<-adjacent(r, N_cells , directions='bishop', pairs=FALSE)

ng_coords <- xyFromCell(r, e)
plot(r)
points(pts, col='red', pch=16)
points(ng_coords, col='black', pch=16)


### Alternative to what Ben suggested!
## Set neighborhood matrix.  Focal cell is 0 while neighbors are 1!!
neigh <- matrix(1L, nrow=5, ncol=5); neigh[3,3] <- 0L
e1<-adjacent(r, N_cells , directions=neigh, pairs=FALSE)
ng_coords1 <- xyFromCell(r, e1)
plot(r);
points(pts, col='red', pch=16);
points(ng_coords, col='black', pch=16);
points(ng_coords1, col='blue', pch=16)
### END

HTH,
Vijay.
On Wed, May 30, 2018 at 4:22 PM, Ben Tupper <btupper at bigelow.org> wrote:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      ASANTOS
    

    

      
      #
    

  


  

      
Thanks so much Vijay and Ben, now it works!!!

 ??????????? Originally, I expect 24 for each 4 point (96 total). 
Apologies, in my example preparation I changed the raster objects 
'land_mask' and 'r' and the idea of my second point [45.12069, 
-88.369745] is to force NA results.

 ??????????? The code below is a solution that I expected:

### Alternative to what Ben suggested!
## Set neighborhood matrix.? Focal cell is 0 while neighbors are 1!!
neigh <- matrix(1L, nrow=5, ncol=5); neigh[3,3] <- 0L
e1<-adjacent(r, N_cells , directions=neigh, pairs=FALSE)
ng_coords1 <- xyFromCell(r, e1)
plot(r);
points(pts, col='red', pch=16);
points(ng_coords, col='black', pch=16);
points(ng_coords1, col='blue', pch=16)
### END