0

I am struggling with the task to set one point in each pixel (of a given raster) in a distance of 3km of my points.

I tried buffering in QGIS and R but cannot find a way to set exactly one point in each pixel. I think it will work with creating a matrix around the points but I couldn´t find a way doing this.

Data: 46 species presence points; Raster Stack with environmental Data Res: 100x100m

enter image description here

raster with presence points (zoomed in)

Crosses: pts; false buffers in purple

presence points: black crosses, buffers(false): purple, backgroundraster: allbands

enter image description here

gridpts

code:

# Points in Matrix around Presence Data 

# #---------------------------------------------------------------------
# Prework

# Packs 
library(spdep)
library(sp)
library(rgdal)
library(raster)
library(RStoolbox)
library(rgdal)

#-------------------------------------------------------------------------------------------------------------------------
# Data 

# Presence Points: 
pts <-readOGR("D://GIDEON_korrigiert_aktuell_repro_EPSG32648_forR.shp")

> pts
class       : SpatialPointsDataFrame 
features    : 46 
extent      : -507313, 812849.9, 761011.2, 2831574  (xmin, xmax, ymin, ymax)
crs         : +proj=utm +zone=48 +datum=WGS84 +units=m +no_defs 
variables   : 2
names       :    Lat,     Lon 
min values  : 10,230, 100,030 
max values  :  7,880,  99,950 

# Raster: 
allbands <- raster("D:/Database/AllReproCliped_for_DOWN/ForestFrek/Stack_8er_Prec_Wet_Temp_Land_Tree_Short_Bare_ForestFrek.tif")

> allbands
class      : RasterLayer 
band       : 1  (of  8  bands)
dimensions : 23913, 17032, 407286216  (nrow, ncol, ncell)
resolution : 106.8938, 106.8938  (x, y)
extent     : -835901.2, 984713.7, 624198, 3180349  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=48 +datum=WGS84 +units=m +no_defs 
source     : D:/Chris/Studium/BA/QGIS_Datein/Database/AllReproCliped_for_DOWN/ForestFrek/Stack_8er_Prec_Wet_Temp_Land_Tree_Short_Bare_ForestFrek.tif 
names      : Stack_8er_Prec_Wet_Temp_Land_Tree_Short_Bare_ForestFrek 
values     : -6.3375, 28.77083  (min, max)

#-------------------------------------------------------------------------------------------------------------------------
# Plot 

plot(allbands)
plot(pts, add= TRUE)

enter image description here

#------------------------------------------------------------------------------------------
# Buffer 

pres_st <- st_as_sf(pts, crs= "bogus")
crs(pres_st)

d = 3000 
buf <- st_buffer(pres_st, dist=d)

gridpts <- st_as_sf(
  as.data.frame(
    mask(allbands, buf),
    xy=TRUE,
    na.rm=TRUE),
  coords=1:2,
  crs=st_crs(allbands))

#------------------------------------------------------------------------------------------
#Plot

plot(gridpts)

enter image description here

plot(allbands)
plot(pts, add=TRUE)
plot(gridpts, add=TRUE)

enter image description here

#------------------------------------------------------------------------------------------


10
  • 1
    I can't quite understand what your expected output is. Your input is 46 points and a raster. Is your output a set of points? 46 of them? Or do you want a point for every pixel in your raster? Across the whole of your raster, or only the non-missing area (land?)?
    – Spacedman
    Apr 27 '21 at 7:44
  • My expected output is a layer with spatial points around my existing presence points. Apr 27 '21 at 8:18
  • Try Raster Sampling in the Toolbox. It samples the raster at the vector locations.
    – wingnut
    Apr 27 '21 at 9:33
  • 1
    Please Edit the question in response to requests for clarification. Code in comments is barely legible to illegible, depending on the language.
    – Vince
    Apr 27 '21 at 13:29
  • 1
    That looks like it has not picked up the X and Y coordinates from the as.data.frame conversion of the raster. Try coords=c("x","y") instead of coords=1:2?
    – Spacedman
    Apr 30 '21 at 13:55
3

I have a dem and some points:

> plot(dem)
> plot(pts, add=TRUE)

enter image description here

I'll make a 1km buffer for my data:

> d = 1000
> buf = st_buffer(pts, dist=d)

Then I can get all the grid points in the buffer by masking the DEM to the buffer and taking all the coordinates of non-NA pixels and making an sf points data frame out of them:

> gridpts = st_as_sf(
      as.data.frame(
        mask(dem,buf),
        xy=TRUE,
        na.rm=TRUE),
      coords=1:2,
      crs=st_crs(dem))

which looks like this:

enter image description here

which is actually nearly 3000 points so close together you can't tell them apart.

4
  • This works perfectly. Thank you very much for that. Unfortunately it does not work with my data. I have updated the post so you can see the problem. When I plot the points by themselves, they are all in line with each other and not near my presence points.I have tried data conversion but could not find a solution. Apr 30 '21 at 8:41
  • If my answer doesnt work then why have you ticked it as an accepted answer? I still can't see the problem even after your edit. What points are "all in a line with each other"? I don't see points in a line on your question.
    – Spacedman
    Apr 30 '21 at 8:49
  • Thanks for noticing that! I have updated the post to make it clearer. Apr 30 '21 at 9:08
  • 1
    Its not clearer unless you include the code you used! From your edits I see you have a SpatialPointsDataFrame which will not work with my code because I'm using the sf package, and functions like st_buffer will not work with that. Show your working in your question alongside the summary outputs of your data objects.
    – Spacedman
    Apr 30 '21 at 9:45
1

Here is an approach with terra

Example data

library(terra)
#terra version 1.2.5
f <- system.file("ex/elev.tif", package="terra")
r <- rast(f)
set.seed(1234)
pts <- spatSample(r, 2, xy=TRUE, values=FALSE, as.points=TRUE, na.rm=TRUE)

Solution:

b <- buffer(pts, 5000)
cnrs <- cells(r, b)
xy <- xyFromCell(r, cnrs[,2])

plot(r)
points(xy)
lines(b, col="red", lwd=2)

enter image description here

2
  • I tried your solution with the exact same data but got the error: Error in .local(x, ...) : not used argument (values = FALSE, as.points = TRUE) when running pts <- .... // terra version: terra_1.1-17 Apr 28 '21 at 15:57
  • You need current terra (version 1.2-5) May 4 '21 at 5:00

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.