3

I'm trying to write code using the Processing R Provider scripting in QGIS. The script is to perform a stratified random sampling scheme on a classified land cover raster image. The land cover codes are 1 (urban), 2 (agriculture), and 3 (water).

The code that I use in R is:

library(sp)
library(rgdal)
library(raster)
landclass <- raster("Land8_sup_class_min_mc.tif")
names(landclass) <- 'stratum'
stratrand <- sampleStratified(landclass, size=10, xp=TRUE, sp=TRUE)
writeOGR(stratrand, dsn = ".", layer = "stratrand", driver = "ESRI Shapefile")

The issue I am having is converting this R code to a syntax that can be used in QGIS. This is one of the many iterations I have tried thus far:

##XRaster = raster
##Size=number 10
##stratrand = output vector
library(sp)
library(raster)
R = raster(XRaster)
values(R) <- 1:ncell(R)
pts = sampleStratified(R, Size, xp=TRUE, sp=TRUE)
df <- data.frame(FID=row.names(pts), row.names=row.names(pts))
pts2 <- SpatialPointsDataFrame(pts, data=df)
stratrand = pts2

I get an error if I run the "names(landclass) <- 'stratum'" portion of the code from R and the values(R) portion of the code runs forever. Any suggestions?

4

The solution:

##landclass=raster
##Size=number 10
##stratrand=output vector
library(sp)
library(rgdal)
library(raster)

names(landclass[[1]]) <- "stratum"
stratrand <- sampleStratified(landclass[[1]],
                              size = Size,
                              xp = TRUE,
                              sp = TRUE)
stratrand

Explanation: QGIS reads any raster as a RasterBrick. If you do raster() on a RasterBrick, you lose the values. Indexing by 1 ([[1]]) takes the first (and the only) layer from a RasterBrick and creates a RasterLayer.

1
  • This was perfect! Thank you for the assistance!
    – Rachel
    Sep 4 '19 at 18:57

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.