Creating a raster layer with a weighted random sample of points (or, my first attempt to create a python script)

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.

The layer with the percentage of cropland IIASA/IFPRI global cropland map; See Fritz S, McCallum I, Schill C, Perger C, See L, et al. (2012) Geo-Wiki: An online platform for improving global land cover. Environmental Modelling & Software 31: 110–123. For the curious, this is east of Lake Victoria around the border between Kenya and Tanzania.

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

The function with as input the weight layer (here the crop mask), the name of the output layer and the minimum and maximum weight. The values of the weight layer should fall within this range (the script does not check this yet)

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.

The resulting binary map with in red the selected raster cells

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:

Percentage of cells with 1 and 0 in each percentage class of the IIASA_crops raster 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.