[R] (no subject)

Sat Oct 7 13:42:56 CEST 2006

On Sat, 7 Oct 2006 Ted.Harding at 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 at nessie.mcc.ac.uk>
> Fax-to-email: +44 (0)870 094 0861
> Date: 07-Oct-06                                       Time: 09:53:43
> ------------------------------ XFMail ------------------------------
> 
> ______________________________________________
> R-help at 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 at nhh.no