9

I want to create a point vector dataset of 10000 points (or larger) within a binary raster, where the points should be constrained to areas where the raster value is 1.

I tried the following steps.

  1. Polygonize raster
  2. QGIS: Vector -> Research Tools -> Random Points

This works fine up to 2000 points but anything above just causes QGIS to crash.

Is there a way to create a vector dataset with a large number of point features constrained by a binary raster (or the polygonized version of it)?

The following tools are at my disposal, ranked from most to least favourable: QGIS, Python, R, ArcGIS

This is what I'm going for, only with 10x the point features.

1k random points

  • How big is your raster, typically? – Spacedman Oct 1 '15 at 9:34
  • The one in the example above is 19200 x 9600. The typical raster is around 10000 x 10000 pixels. – Kersten Oct 1 '15 at 10:51
  • Okay, the more RAM your machine has the better. I dare not test on a 10,000x10,000 raster on my little PC here, although you could always divide the raster, sample in parts, and join... – Spacedman Oct 1 '15 at 10:57
  • why polygonize the raster? do you mind this answer is useful for you? gis.stackexchange.com/questions/22601/… – Luigi Pirelli Oct 2 '15 at 9:10
  • Because then I can use the "Random Points in Polygon" function, whereas QGIS doesn't have a "Random Points inside specific values of a Raster" function. – Kersten Oct 2 '15 at 10:58
7

Here's a way in R:

Make a test raster, 20x30 cells, make 1/10 of the cells set to 1, plot:

> require(raster)
> m = raster(nrow=20, ncol=30)
> m[] = as.numeric(runif(20*30)>.9)
> plot(m)

For an existing raster in a file, for example a geoTIFF, you can just do:

> m = raster("mydata.tif")

Now get a matrix of the xy coordinates of the 1 cells, plot those points, and we see we have cell centres:

> ones = xyFromCell(m,1:prod(dim(m)))[getValues(m)==1,]
> head(ones)
       x    y
[1,] -42 85.5
[2,] 102 85.5
[3,] 162 85.5
[4,]  42 76.5
[5,] -54 67.5
[6,]  30 67.5
> points(ones[,1],ones[,2])

Step 1. Generate 1000 (xo,yo) pairs that are centred on 0 in a box the size of a single cell. Note use of res to get the cell size:

> pts = data.frame(xo=runif(1000,-.5,.5)*res(m)[1], yo=runif(1000,-.5,.5)*res(m)[2])

Step 2. Work out which cell each of the above points is going into by randomly sampling 1000 values from 1 to the number of 1 cells:

> pts$cell = sample(nrow(ones), 1000, replace=TRUE)

Finally compute the coordinate by adding the cell centre to the offset. Plot to check:

> pts$x = ones[pts$cell,1]+pts$xo
> pts$y = ones[pts$cell,2]+pts$yo
> plot(m)
> points(pts$x, pts$y)

Here's 10,000 points (replace the 1000 above with 10000), plotted with pch=".":

points in ones

Pretty much instantaneous for 10,000 points on a 200x300 raster with half the points as ones. Will increase in time linearly with how many ones in the raster, I think.

To save as a shapefile, convert to a SpatialPoints object, give it the right coordinate system reference (the same as your raster) and save:

> coordinates(pts)=~x+y
> proj4string(pts)=CRS("+init=epsg:4326") # WGS84 lat-long here
> shapefile(pts,"/tmp/pts.shp")

That will create a shapefile that includes the cell number and offsets as attributes.

  • This is looking very promising. My R has gotten a little rusty: How would I be able to export the points to a vector format (Shapefile, geojson, gml,...whatever really) - I need to save the locations of the sample points for later use. – Kersten Oct 1 '15 at 10:49
  • Edits show how to read a raster and convert pts to shapefile... – Spacedman Oct 1 '15 at 10:56
3

Whenever I work with large datasets, I like to run tools/commands outside of QGIS such as from a standalone script or from OSGeo4W Shell. Not so much because QGIS crashes (even if it says "Not responding", it's probably still processing the data which you can check from the Task Manager), but because more CPU resources such as RAM are available to process the data. QGIS itself consumes a fair chunk of memory to run.

Anyway, to run a tool outside QGIS (you would need to have installed QGIS via the OSGeo4W installer), follow the first 2 steps as described by @gcarrillo in this post: Problem with import qgis.core when writing a stand-alone PyQGIS script (I suggest to download and use his .bat file).

Once the PATHS are set, type python into the command line. For convenience, copy the following code into a text editor such as Notepad, edit the parameters such as the pathname of your shapefile etc. and then paste the whole thing into the command line by Right-click > Paste:

import os, sys
from qgis.core import *
from qgis.gui import *
from PyQt4.QtGui import *

from os.path import expanduser
home = expanduser("~")

QgsApplication( [], False, home + "/AppData/Local/Temp" )

QgsApplication.setPrefixPath("C://OSGeo4W64//apps//qgis", True)
QgsApplication.initQgis()
app = QApplication([])

sys.path.append(home + '/.qgis2/python/plugins')
from processing.core.Processing import Processing
Processing.initialize()
from processing.tools import *

shape = home + "/Desktop/Polygon.shp"
result = home + "/Desktop/Point.shp"
general.runalg("qgis:randompointsinlayerbounds", shape, 10000, 0, result)

Using the script, I ran the Random points in layer bounds tool for a fairly large shapefile and it took under 20 seconds to generate 10k points. Running it inside QGIS took almost 2 minutes so atleast for me, there's a significant difference.

  • 1
    Excellent alternative, +1. Just tested this for my application and while it is a bit slower than the R approach it creates the desired results. – Kersten Oct 2 '15 at 12:15
  • @Kersten - Awesome, glad it works :) – Joseph Oct 2 '15 at 12:19
1

You can also use GRASS GIS directly for this job - Stratified random sampling: Random sampling from vector map with spatial constraints:

https://grass.osgeo.org/grass72/manuals/v.random.html#stratified-random-sampling:-random-sampling-from-vector-map-with-spatial-constraints

Additionally, random sampling from vector map by attribute and a few other methods are implemented in the command.

Note: The v.random version exposed in QGIS through processing does not reflect the full functionality but just a simplified view.

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.