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.


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
do.call("rbind", spdf) -> spdf

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, as.data.frame(wrld_simpl), 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.

  • 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){
  # 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,nrow(x))
  x[c('geoX','geoY')] = pts
  # returning the completed chunk

I use the term "geocode" extremely loosely here. What this does, is you pass in all the data for a particular polygon, as well as the polygon set. (would be faster to pass just the polygon, but i didn't optimize to that extent). It computes random XY coordinates, then attaches those XY coordinates to the data and returns it to the main program. "data" (in the complete sense) 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.

  • I'm trying to understand this and get it to work. Should the first line of the function be: cur.poly <- polygons[polygons$polyID == x[1,]$polyID,] ? (ie. missing =) Also, where does data$polyID come from in data = data[(data$polyID %in% shapeFile$polyID),] ? I have no clue on that one. – Simon Apr 28 '17 at 21:13
  • This code was obfuscated from working code... i think I missed a couple things in doing the obfuscation. Yes, that should be a ==, as it's building a boolean vector. In truth, this is rather inefficient code and i'd do it a slightly different way were i doing it now. – Faydey Apr 29 '17 at 21:57
  • data$polyID refers to the polygon ID of the polygon you want the XY coords generated in. data refers to the data you want to scatter. in my case, i was doing this for HMDA data. polyID refers to census tract, identified by 11 digit FIPS. Sometimes, the reported FIPS for hmda data is not valid, and i would get show-stopping errors (function would crash) when that was the case. that step is to ensure that only data with valid FIPS gets through. – Faydey Apr 29 '17 at 22:05
  • Ok, so if I want random points in all polygons then I just order and split the shapeFile$polyID (in this example) and ignore thedata dataframe? – Simon Apr 30 '17 at 8:14


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?

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