[R] random point in a circle centred in a geographical posit

[(Ted Harding)]

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.

> 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 at nessie.mcc.ac.uk>
Fax-to-email: +44 (0)870 094 0861
Date: 07-Oct-06                                       Time: 09:53:43
------------------------------ XFMail ------------------------------