4

I have a projected raster shown in the image below (left). I want to generate random points inside the raster like the middle image. However, the result I am getting is the same as the image on the right side. How would I get the random points inside the extent of the projected raster?

enter image description here

library(rgdal)
library(raster)

a <- brick("E:/raster.tif")

b <- as(extent(a), "SpatialPolygons")

c <- spsample(b, n = 20, type = "random")

par(mar=c(1,1,1,1))
plot(c)
3

Suppose a Landsat band:

r <- raster('~/path/to/LC80160442016191LGN00_cfmask_conf.tif')

plot(r)

enter image description here

I recommend to compute NA cells (to determine raster's edges) and convert them to lines (because polygonize is a very slow process in R). After this, convert them to polygon and use it to sample points:

lin <- rasterToContour(is.na(r))

library(sf) # to call st_* functions

pol <- as(st_union(st_polygonize(st_as_sf(lin))), 'Spatial') # st_union to dissolve geometries

pts <- spsample(pol[1,], 2000, type = 'random')

plot(r)
plot(pts, add = T, col = 'red')

enter image description here

  • Wow. It worked! Thank you very much @aldo_tapia. Now, how would I grab the pixel values of those 2000 random points so I could look at them separately (maybe through spreadsheets)? – GISnew Aug 3 '18 at 15:57
  • @GISnew with rate you mean sort them? Because with data.frame(pts@coords, pts@data) you can extract coordinates and data.frame related (after extract values) – aldo_tapia Aug 3 '18 at 16:07
  • From the initial code you provided, I actually exported 'pts' as .csv and saw that it already contains the x and y coordinates. However, I would want a third column with the value of the pixels of the Landsat band from those points. – GISnew Aug 3 '18 at 17:42
  • @GISnew Use extract() function to intersect points with raster. You had not defined this in your question – aldo_tapia Aug 3 '18 at 18:38
  • Thanks @aldo_tapia. Got it. I created a data.frame then this: pts2 <- cbind(extract(r, xy, df = T),xy). – GISnew Aug 3 '18 at 20:42
3

Rasters are always rectangles. Your shape on the far left is padded with no-data cells, so the extent object will be a rectangle that contains the visible data and the NA padding.

To sample randomly from data cells only, you could try something like

# which cells are not NA? These ones:
notna <- which(!is.na(values(raster)))

# grab 20 cell index numbers at random
samp <- sample(notna, 20, replace = FALSE)

# and their values
sampdata <- raster[samp]

# and their location coordinates
samplocs <- xyFromCell(raster, samp)

# convert to a data frame
samp <- as.data.frame(cbind(samplocs, samp, sampdata))
names(samp) <- c('x', 'y', 'index', 'value')

# and optionally, a spatial object, in sf or sp
library(sp)
samp_sp <- samp
coordinates(samp_sp) <- ~x + y
crs(samp_sp) <- CRS('+init=epsg:4326')

library(sf)
samp_sf <- st_as_sf(as.data.frame(samp), coords = c('x', 'y'), crs = 4326)

This does preclude sampling more than once from any given cell, and all sample points will be in the exact middle of cells, which may or may not matter to you.

  • Thanks @obrl_soil. I tried the code but when I run sampdata <- raster[samp] Error in raster[samp] : object of type 'closure' is not subsettable – GISnew Aug 2 '18 at 17:50
  • what do your samp and raster objects look like when you print them to the console? (oh wait you might have to specify base::sample() depending on what other packages are loaded - functions sometimes have conflicting names) – obrl_soil Aug 2 '18 at 21:42
  • Also, if your raster is actually multiband (I see now you've used brick()), you will need to change your code to which(!is.na(values(raster[[1]]))). This will assume that all of your raster bands have the same no-data areas. – obrl_soil Aug 2 '18 at 22:03
  • thank you for the help @obrl_soil. I finally resolved it. – GISnew Aug 3 '18 at 20:45

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.