Sign up ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I have an area divided into polygons, and a dataset of observations which has a common field with the polygon layer.

I now would like to generate a random point for each dataset entry within the according polygon.

R was already a great help in digesting the dataset and even to produce shapefiles with summed values of the dataset, so I dare to post this here in the gis.stackexchange. It would be nice to have a solution in R, cause i have to create multiple subsets which can so easily be done there by a loop, so i do no have to filter the data manually in QGIS all the time.

Other solution are welcome as well. enter image description here

Polygonlayer: 39 Polygons


Datalist: 1360 Data entries

Polyid, value1, value2, value3
A17, ..., ..., ...
O6, ..., ...
A3, ...
A4, ...
O6, ...
A8, ...
A17, ...

Goal: -> Having a point-shape with 1360 random points, which reside in the according polygon and have all attributes of the datalist.

share|improve this question

4 Answers 4

up vote 1 down vote accepted

See my example and answer

### Preparing the SpatialPointsDataFrame
spdf <- matrix(as.numeric(NA), nlevels(Poly$MatchID), 1)
spdf <- as.list(spdf)

### Sample the coordinate, match it with data in spdf. It create a list fore each factor of the MatchID
### sample(spsample()) fix the size of the sample
for (i in seq(Poly$MatchID))
    spdf[i] <- SpatialPointsDataFrame(
               sample(spsample(Poly[order(Poly$MatchID)==i,], n = 100, "stratified"),table(df$MatchID)[[i]]),  ### table(df$MatchID)[[i]] is the size of the sample and match the sum of factors in df 
               df[df$MatchID==dimnames(table(df$MatchID))[[1]][i],], ##  dimnames(table(df$MatchID))[[1]][i] ### match the value of the selected "factor" to select the rows of the data

## Merging together the list to make on SpatialDataFrame"rbind", spdf) -> spdf
share|improve this answer

You can use spsample for this. You need to set the "iter" argument depending on how your polygons are shaped, but this code generates a coordinate from each polygon in the wrld_simpl data set. Using iter=10 found a coordinate for each polygon for these data (but 4 failed a few times).

s <- matrix(as.numeric(NA), nrow(wrld_simpl), 2)
for (i in seq_len(nrow(s))) s[i,] <- as.vector(coordinates(spsample(wrld_simpl[i,], type = "random", n = 1, iter = 10)))

Balancing the right ways of sampling in order to ensure always an intersecting point is always found may take some care, see the help page details for ?spsample. Note this is sampling in longlat so it's not exactly sensible, but that will depend on your data, not the function.

This is treating each row in the SpatialPolygonsDataFrame as an entity to sample from, so only one point per "multipolygon". There are methods for Polygon, Polygons and SpatialPolygons so it will work similarly for those with slightly different code depending on what you have.

Note that coordinates(x) will give you centroid of each polygon, but that's obviously not guaranteed to land in the shape itself. Also, spsample assumes that any holes are sensibly defined.

With that matrix each row is a one-to-one match to rows in the polygon layer, so

res <- SpatialPointsDataFrame(s,, proj4string = CRS(proj4string(wrld_simpl)))

finishes the job, but it depends on what your objects are and how you deal with them. Record matching and transfer is pretty basic in R, but not exactly on-topic here.

share|improve this answer
sorry for answering earlier. had to dig deeper into R (great success so far). Unfortunately, i can't achieve what I want to do with your proposal. I only get as many points as I have polygons, but i need as many points as i have data rows in my list. My idea was to loop through the list, pick the Polygon-ID of the first row, select the according polygon in my shape with the same ID, create a random point, attach the rest of the list data to it ... next row ... – Bernd V. Oct 24 '13 at 14:49
Right, I missed that part. You need to replace wrld_simpl[i, ] with wrld_simpl[s[i,"poly-id", ] or similar. Reproducible examples help, I'll try to make one for you if I get around to it. – mdsumner Oct 24 '13 at 21:10

I recently needed to do something similar. I had around 70,000 polygons, and around ~40 million points to place. With that in mind, a few things about the particulars of my answer:

  1. it needs to be fast.
  2. nothing else is relevant to the placement of the point other than the polygon ID.
  3. none of the polygons I am working with are nested inside any signle polygon. (More accurately, none of the polygons i'm working with have holes in them.

Given these, I was able to make use of the csr() function from splancs. All i need to do is pass a bounding box (the polygon in question, using vertex coordinates list), and the number of points to generate, and it will have at.

library(snowfall) # for performance
library(rgdal)    # for the shapefile storage

So I'm making use of the snowfall package for parallelization. I use 4 threads. I could have used all 8, but the size of the dataset is so large I would have ran out of memory.

geocode = function(x,polygons){
  # some housekeeping = x[1,]$polyID
  n = nrow(x)
  # creating a target to pass data into
  x[c('geoX','geoY')] = NA
  # selecting the polygon, and getting its coordinates
  cur.poly = polygons[polygons$polyID = x[1,]$polyID,]
  boundingBox = cur.poly@polygons[[1]]@Polygons[[1]]@coords
  # Generating the points, passing them back into the target.
  pts = csr(boundingBox,n)
  x[c('geoX','geoY')] = pts
  # returning the completed chunk

I use the term "geocode" extremely loosely here. "data" refers to the 40 million observation table that I want to attach XY coordinates to. I subset the data table because there were a few observations with malformed polygon ID's (polygon ID's which didn't appear in the polygon shapefile). Leaving them in caused an error, so better to remove.

# Using readOGR from rgdal package
shapeFile = readOGR(dsn = shapefilePath, layer = shapefileName)
# subsetting ensure that the data has a valid polygon to go into
data = data[(data$polyID %in% shapeFile$polyID),]
# splitting the data into a list 
data = data[order(data$polyID),]
chunks = split(data,data$polyID)

So the way snowfall works is kind of interesting. It spawns child R processes, that are blank slates. You need to explicitly pass to them beforehand whatever you intend to work with, that doesn't get passed by "apply" (basically, everything but x). So here, that means passing the function definition, along with the shapefile. As well as importing the packages we'll be using.

# initializing the cluster
sfInit(parallel = TRUE, cpus=4)
# some housekeeping

All this reduces to 1 line of code. This is, however, insanely fast for this purpose. I had done this batch scattering in polygon before using arcGIS, and it took the better part of the weekend to do all 40 million observations. Subsetting by polygon then doing this scatter had the whole deal done in ~15 minutes.

# Where the magic happens
output = sfLapply(chunks,function(x){geocode(x,shapeFile)})

At the end we are left with a list of identically schema-ed data frames. I actually still don't know a clever way to stack these, but I didn't need it in my workflow. (I pushed them back to an SQLite database, appending onto the end of the table).

I realize this is long after this question was answered, but hopefully others find and like this.

share|improve this answer


2nd attempt:

How about Vector -> "Random points" selecting "Poyid" field on "value from input field" and then "Intersect" to transfer the attributes from one layer to the other?

share|improve this answer
well, i tried to explain it, but obviously didn't reach the goal. Generating 1360 points "somewhere" on the layer does not do the trick. I need e.g. 5 points in A1, 250 points in 017 etc. I search for a way to take the data entry, tell some process to extract the polyid, generate a random point within this polygon, merge the rest of the attributes to the point -> next data entry etc. – Bernd V. Oct 13 '13 at 16:09
I don't think you can do all that in one step. But see my second answer. – Filipe Dias Oct 13 '13 at 16:23

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.