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.

  • Do you need exactly that percentage of pixels selected within each polygon or do you want to select each pixel independently with the given probability? The latter is much easier to do--and corresponds to many common random sampling procedures--but will typically select a slightly different number than actually targeted. – whuber Jan 16 '13 at 15:07
  • 1
    @whuber I'm looking for the exact solution (as far as possible with integer number of pixels). – sebschub Jan 16 '13 at 15:14
  • Maybe better asked on the user list: grass.osgeo.org/support/mailing-lists since way more people are reading there. – markusN Jan 16 '13 at 23:21

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:

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="|")

## 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.

| improve this answer | |

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.