From: Roger Bivand <Roger.Bivand_at_nhh.no>

Date: Sat 07 Oct 2006 - 11:42:56 GMT

Date: Sat 07 Oct 2006 - 11:42:56 GMT

On Sat, 7 Oct 2006 Ted.Harding@nessie.mcc.ac.uk wrote:

> Hi Albert

*>
**> On 07-Oct-06 Albert Picado wrote:
**> > Dear List members
**> >
**> > I am trying to find a way to generate a random point in a circle
**> > centred in a geographical location.
**> >
**> > So far I have used the following formula (see code below):
**> > random_x = original_x + radius*cos(angle)
**> > random_y = original_y + radius*sin(angle)
**> > where radius is a random number between 0 and the radius of the
**> > circle and angle is between 0 and 360 degrees
**> >
**> > The code bellow works fine however some people argue that this method
**> > will produce uneven locations--I really do not understand why.
**>
**> Think about 5000 points chosen by the method you describe, and
**> divide the circle into 10 circular regions at r = 0.1, 0.2, ... 1.0
**> at equal increments of r.
**>
**> Since you have chosen r uniformly distributed (in your description
**> above and in your code below) you will have (about) 500 points in
**> each region. So that is 500 between r=0 and r=0.1, in an area
**> pi*(0.1^2); and then 500 between r=0.1 and r=0.2, in an area
**> pi*(0.2^2) - pi*(0.1^2), and so on.
**>
**> The second area is pi*(0.04 - 0.01) = 0.03*pi, which is 3 times
**> the first area, 0.01*pi. And so on. So the average density of
**> points in the first is 3 times the average density of points in
**> the second. And so on.
**>
**> Therefore your method will give points whose density per unit
**> area is more concentrated the nearer you are to the centre of
**> the circle.
**>
**> You can check this visually with the following code:
**>
**> r <- runif(5000)
**> t <- 2*pi*runif(5000)
**> x <- r*cos(t) ; y <- r*sin(t)
**> plot(x,y,pch="+",col="blue",xlim=c(-1,1),ylim=c(-1,1))
**>
**> The above argument should suggest the correct way to proceed.
**> You need equal increments in "probability that radius < r"
**> to correspond to equal increments in the area of a circle with
**> radius r. This means that the "probability that radius < r"
**> should be proportional to r^2. Hence make the random radius r
**> to be such that r^2 is uniformly distributed.
**>
**> Now try the following modification of the above code:
**>
**> r <- sqrt(runif(5000))
**> t <- 2*pi*runif(5000)
**> x <- r*cos(t) ; y <- r*sin(t)
**> plot(x,y,pch="+",col="blue",xlim=c(-1,1),ylim=c(-1,1))
**>
**> and you will see that (to within random scatter) the result is
**> uniformly distributed over the circle. (By the way, don't use
**> degrees from 0 to 360 in sin() and cos() -- these use radians
**> from 0 to 2*pi).
**>
**> > Another option will be to use the “rejection method: generate
**> > points inside a box and reject points that fall outside the circle.
**> > However I think this is much complicated--or at least I don't know
**> > how to programme it.
**>
**> It is more complicated. Uniform x and y are easy, but then you have
**> to check that x^2 + y^2 <= 1 and reject it if not; and also keep
**> a count of the number of points you have accepted and check whether
**> this has attained the number of points you need, so that you can
**> continue sampling until you have enough points (which you cannot
**> predict at the start).
**>
**> It is also (though in this case not very) inefficient. A unit
**> circle has radius pi = 3.14... and the enclosing square has
**> area 4, so you will reject nearly 25% of your points.
**>
*

Since these were geographical points, I'd risk using the sp package for this, but first a little function to make a polygon-circle based on Gabor Csardi's answer to a similar question in April:

circle <- function(x, y, r, length=100, zapsmall=TRUE) {
t <- seq(0, 2*pi, length=length)

coords <- t(rbind(x+sin(t)*r, y+cos(t)*r))
if (zapsmall) coords <- zapsmall(coords)
coords

}

Then:

library(sp)

crds <- circle(0, 0, 1)

Poly_0.0.1 <- Polygon(crds)

set.seed(061007)

pts <- sp:::sample.Polygon(Poly_0.0.1, n=100, type="random")
plot(crds, type="l", asp=1)

dim(coordinates(pts))

points(coordinates(pts))

which uses rejection. The distances follow as:

spDistsN1(coordinates(pts), c(0,0))

Roger

> > So,

*> > 1. Is there any package that generates random points in a circle? (I
**> > couldn’t find any)
**>
**> Use the above code as a basis. For n points in a circle of radius
**> 1000, multiply sqrt(runif(n)) by 1000.
**>
**> > 2. Do you think the code bellow produces uneven locations?
**>
**> Yes (see above)
**>
**> > 3. Any suggestions to programme the "rejection method"? Maybe someone
**> > has already used it...
**>
**> I wouldn't bother ... ! But if you are really interested in how
**> it can be done, please write again.
**>
**> > Cheers
**> >
**> > albert
**> >
**> >#Code random location in a circle
**> ># Generate UTM coordinates (original location)
**> >> x_utm<-c(289800)
**> >> y_utm<-c(4603900)
**> >> coord<-as.data.frame(cbind(x_utm,y_utm))
**> ># Generate a random radius with a maximum distance of 1000 meters from
**> ># the original location.
**> >>radius<-runif(1,0,1000)
**> ># Generate a random angle
**> >>angle<-runif(1,0,360)
**> >#Generate a random x coordinate
**> >>x_rdm<-coord$x_utm[1]+radius*cos(angle)
**> >#Generate a random y coordinate
**> >>y_rdm<-coord$y_utm[1]+radius*sin(angle)
**> >>result<-as.data.frame(cbind(x_rdm,y_rdm, x_utm, y_utm))
**> ># We can calculate the distance between the original and the random
**> ># location
**> >> result$dist_m<-sqrt(((coord$x_utm-result$x_rdm)^2+
**> > (coord$y_utm-result$y_rdm)^2))
**> >>result
**>
**> Best wishes,
**> Ted.
**>
**> --------------------------------------------------------------------
**> E-Mail: (Ted Harding) <Ted.Harding@nessie.mcc.ac.uk>
**> Fax-to-email: +44 (0)870 094 0861
**> Date: 07-Oct-06 Time: 09:53:43
**> ------------------------------ XFMail ------------------------------
**>
**> ______________________________________________
**> R-help@stat.math.ethz.ch mailing list
**> https://stat.ethz.ch/mailman/listinfo/r-help
**> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
**> and provide commented, minimal, self-contained, reproducible code.
**>
*

-- Roger Bivand Economic Geography Section, Department of Economics, Norwegian School of Economics and Business Administration, Helleveien 30, N-5045 Bergen, Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43 e-mail: Roger.Bivand@nhh.no ______________________________________________ R-help@stat.math.ethz.ch mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.Received on Sat Oct 07 21:49:28 2006

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.1.8, at Sat 07 Oct 2006 - 13:30:10 GMT.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help.
Please read the posting
guide before posting to the list.
*