I needed to create a raster map layer with a weighted random sample of all raster cells, using the percentage of crop land as weight. I couldn’t find a function to create such a weighted sample, so I decided to create a script to do this for me.
Update: I created an addon based on the script presented in this post. It is available through the GRASS addons svn repository.
The script (here) basically creates a new layer with random numbers between a minimum and maximum value. Next it checks whether the random number is smaller than the weight of the corresponding cell in an user defined raster layer (the weights layer). If so, a number 1 is assigned, otherwise a 0.
Very basic, only a few commands.. in other words, ideal to finally get out of my comfort zone (R or shell script, sort of) and try to do this with a python script. And, using some existing scripts as example and information from this wiki page, it turned out to be easier than I thought.
Below some screenshots where I use a cropland layer (percentage of cropland, see image above) as weight to create a new layer with sample points (the script is at the end of this post).
And below the resulting layer. In red the raster cells that were assigned a value 1, in white the cells with value 0. If you need to select fewer points, you can use the layer below as a mask and run r.random to get a subset of these raster points.
To check whether the result really represents a weighted sample, one can for example run
r.report -h -n --overwrite map=IIASA_crops,testmap units=p to get the percentage of selected points for each unique value in the weight layer:
It is a very simple script of course, but it shows nonetheless that with a obviously very rudimentary/ hardly existing knowledge about Python, one can still create simple scripts in GRASS GIS ;-). Next step would be for the script to check whether the minimum and maximum values of the weight layer fall within the range defined by the start and end values provided by the user (any suggestion?).
There also should be the option (or default) that the script uses the minimum and maximum value of the weight layer rather then needing the user to supply a minimum and maximum. Thinks I know (I think) how to do in a shell script, but not yet in Python.
[Update: The addon now compares the user provided minimum and maximum weight with those of the weight raster, and gives a warning if the user defined minimum > minimum raster value or if the user defined maximum < maximum raster value. It is now furthermore possible to set the number of sample points (as absolute number or percentage)]
And how to get a progress indicator? There is something about that on the wiki but I do not really know how to use that. Anyway, this little exercise was enough to show me that it might be worth it spending a bit of time on scripting in Python.. now if I only could find that time. You can check out the code here.