[R] Re : Draw a circle with a given radius in an existing map

[Arien Lam]

A function that works for me, assuming the Earth is close 
enough to a sphere (which may or may not be true in your 
application), follows below.

Hope this helps, Arien

# computes the place where you end up, if you travel a 
certain distance along a great circle,
# which is uniquely defined by a point (your starting point)
# and an angle with the meridian at that point (your direction).
# the travelvector is actually a dataframe with at least the 
colums
# magnitude and direction.
# n.b. earth radius is the "ellipsoidal quadratic mean 
radius" of the earht, in m.

vectordestination <- function(lonlatpoint, travelvector) {
     Rearth <- 6372795
     Dd <- travelvector$magnitude / Rearth
     Cc <- travelvector$direction

     if (class(lonlatpoint) == "SpatialPoints") {
         lata <- coordinates(lonlatpoint)[1,2] * (pi/180)
         lona <- coordinates(lonlatpoint)[1,1] * (pi/180)
     }
     else {
         lata <- lonlatpoint[2] * (pi/180)
         lona <- lonlatpoint[1] * (pi/180)
     }
     latb <- asin(cos(Cc) * cos(lata) * sin(Dd) + sin(lata) 
* cos(Dd))
     dlon <- atan2(cos(Dd) - sin(lata) * sin(latb), sin(Cc) 
* sin(Dd) * cos(lata))
     lonb <- lona - dlon + pi/2

     lonb[lonb >  pi] <- lonb[lonb >  pi] - 2 * pi
     lonb[lonb < -pi] <- lonb[lonb < -pi] + 2 * pi

     latb <- latb * (180 / pi)
     lonb <- lonb * (180 / pi)

     cbind(longitude = lonb, latitude = latb)
}


# You can use this function to draw a circle with radius of 
ten kilometer around a point:

radius <- 10000 # in meter
NewHaven <- c(-72.9,41.3) #c(lon,lat)
circlevector <- as.data.frame(cbind(direction = seq(0, 2*pi, 
by=2*pi/100), magnitude = radius))

mycircle <- vectordestination(NewHaven, circlevector)

plot(mycircle,type="l")





Roger Bivand schreef:
> On Tue, 7 Nov 2006, justin bem wrote:
> 
>> Try this :
>>     >it<-seq(0,2*pi, l=100)
>>     >xt<-r*cos(it)
>>     >yt<-r*sin(it)
>>     >points(xt,yt,type="l",col="blue")
>>
>> a circle of radium r is define by
>>    xt=r*cos(t)
>>    yt=r*sin(t) 
> 
> Isn't this suggestion on the plane, when the question was about finding
> the coordinates on the surface of the sphere (globe) in degrees of
> longitude and latitute that are x miles from the centre?
> 
> If the area is not large, then projecting the centre point to a suitable 
> planar projection, making the circle on the plane as above, and inverse 
> projecting back to geographical coordinates should work (function 
> project() in package rgdal). If the radius is in thousands of miles, the 
> projection distortion would be considerable, though.
> 
>> Justin BEM
>> Elève Ingénieur Statisticien Economiste
>> BP 294 Yaoundé.
>> Tél (00237)9597295.
>>
>>
>>
>> ----- Message d'origine ----
>> De : Xiaomei Ma <xiaomei.ma at yale.edu>
>> À : R-help at stat.math.ethz.ch
>> Envoyé le : Mardi, 7 Novembre 2006, 6h11mn 27s
>> Objet : [R] Draw a circle with a given radius in an existing map
>>
>>
>> I have drawn a map in which the X and Y axes are latitude and 
>> longitude. Now I need to draw one circle on the map - the center is a 
>> point with specific latitude and longitude, but the challenge is that 
>> the radius is in miles. Is there a way to do this? I'd very much 
>> appreciate your response.
>>
>> XM
>>
>> ______________________________________________
>> 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.
>>
>>
>> 	
>>
>> 	
>> 		
>> ___________________________________________________________________________ 
>> Découvrez une nouvelle façon d'obtenir des réponses à toutes vos questions ! 
>> Profitez des connaissances, des opinions et des expériences des internaut
>>
>> 	[[alternative HTML version deleted]]
>>
>>
> 

-- 
drs. H.A. (Arien) Lam (Ph.D. student)
Department of Physical Geography
Faculty of Geosciences
Utrecht University, The Netherlands

E-mail:     a.lam at geo.uu.nl
Web:        http://www.geo.uu.nl/staff/a.lam

Telephone:  +31(0)30-253 2988 (direct), 2749 (secretary)
Fax:        +31(0)30-2531145

Visiting address: Room Zon-121, Heidelberglaan 2, Utrecht
Postal address: P.O.Box 80.115, 3508 TC Utrecht