Set certain amount of pixels to 1 inside of a polygon (GRASS)

Question

I have a rather coarse vector map with polygons and want to calculate a (sufficiently) finer raster map defined by spatial extension and resolution. Every polygon of the vector map has a value associated that describes the percentage of pixels of the finer map inside the polygon that should be randomly set to 1. The fraction of pixels set to 1 should be as exact as possible.

Example: Polygon 7 is associated with a value of 20%. 50 pixels of the finer map would be (mainly) situated inside Polygon 7. Thus, 50*0.2=10 out of the 50 should be set to 1.

As a possible approach in GRASS, I found v.to.rast with use=cat to create a raster map from the polygons where the value indicates in which polygon they lie. Now, I would need to select randomly a percentage of the pixels with the same value.

Answer 1

This can be solved in conjunction with R. First, as assumed in the question, v.to.rast with use=cat is used to find the pixels of the finer raster map that lie in each polygon. Subsequently, R is used to set the correct amount of pixels to 1.

The following scripts have to be adopted to the specific needs. In the driving GRASS script, irr is the vector map of polygons. Its connected database has the column cat, which identifies each polygon uniquely, and the column IRRI_PCT that gives the percentage of fine pixels that should be set to 1.

Main GRASS script:

# target region and resolution
g.region n=2.4 s=-2.6 w=-2.6 e=3.8 nsres=0.01 ewres=0.01
# find pixels of target map which lie in each polygon of irr
v.to.rast --o in=irr out=irr_rast use=cat

# run R
R CMD BATCH apply_irr.R

R script:

library(spgrass6)
rast <- readRAST6("irr_rast")        # raster with values to use with
                                     # one irrigation fraction
irrrast <- rast                      # target map
irrrast$irr_rast <- 0.               # initialize it

## read table of irrigation fraction
temp_irr <- execGRASS("db.select", parameters=list(table="irr"), intern=TRUE)
temp_con <- textConnection(temp_irr)
irrtab <- read.table(temp_con, header=TRUE, sep="|")
close(temp_con)

## set irrtab$IRRI_PCT % pixels to 1
for (nr in 1:nrow(irrtab)) {
  pixels <- which(rast$irr_rast==irrtab$cat[nr])
  irrrast$irr_rast[sample(pixels, floor(length(pixels)*irrtab$IRRI_PCT[nr]/100.)  )] <- 1.
}

## write results
writeRAST6(irrrast, "irr_rast", overwrite=TRUE)

The final raster map is available in GRASS with the name irr_rast.