3

I would like to use R and leaflet to draw a circle (a real circle on the surface of the earth) on a mercator projection.

In the following code, I have the center of the circle (with lng=0 and lat=60), and I would like to have a parametric equation to draw the circle. I started with the equation below, and I can draw it on a leaflet map:

centre=c(0,60)
t=seq(0,2*pi,0.005)
r=9
xx=cos(t)*r+centre[1]
yy=sin(t)*r+centre[2]

Then I can draw the form on a leaflet map: 

leaflet() %>% 
  addTiles() %>% 
  addPolygons(
    lng=xx, lat=yy,fillColor = "transparent"
  )%>% 
  addCircleMarkers(lng=centre[1], lat=centre[2],radius=1,color="red")

The code above gives this kind of formes on the earth: http://www.mylovedone.com/image/solstice/sum12/tissotFig3b.jpg

I did't make the image, but I think they are the same. And what I want to draw is this: https://i.stack.imgur.com/IiN3c.jpg

How should I modify the equation of the circle so that I can obtain a real circle.

Is this a Tissot's indicatrix? It seems that a Tissot's indicatrix is to render a circle on the sphere on a map. However, the Tissot's indicatrix on a mercator projection is a circle! I can't understand how it is possible. Since if you draw a circle on the north pole, there is no way it is a circle on mercator projection.

EDIT: I found this web site http://www.movable-type.co.uk/scripts/latlong.html it shows the equation of destinations, givens one point, the direction and the distance. So it is then possible to draw the cercle!

John Smith
  • 389
  • 3
  • 10
  • I found this: http://adsabs.harvard.edu/full/1932JRASC..26..161G – John Smith Jun 01 '16 at 08:43
  • (1) Tissot's indicatrix is an infinitesimal circle that has been magnified to be visible on the map. See http://gis.stackexchange.com/questions/5068 for details. (2) Please indicate whether you want the formula in terms of (lat, lon) or in Mercator coordinates: your post and your comments give contradictory indications. – whuber Jun 02 '16 at 19:06

1 Answers1

1

One way is to draw in Orthographic and reproject. It's not very generalizable though, and I'm unsure how accurate this is at 1000skm scales. (Looks good for this example compared to Manifold's geog-circle). 

## centre in longlat
llcentre=c(0,60)

library(rgdal)

ocentre <- c(0, 0)
t=seq(0,2*pi,0.005)
r=2500000  ## this must be in the units of the map projection, i.e. metres
xx=cos(t)*r+ocentre[1]
yy=sin(t)*r+ocentre[2]


sp <- SpatialPoints(cbind(xx, yy), proj4string = CRS(sprintf("+proj=ortho +ellps=WGS84 +lon_0=%f +lat_0=%f", llcentre[1], llcentre[2])))
library(maptools)
data(wrld_simpl)
map <- project(coordinates(as(as(wrld_simpl, "SpatialLines"), "SpatialPoints")), "+proj=merc +ellps=WGS84")

plot(map, asp = 1)
plot(spTransform(sp, "+proj=merc +ellps=WGS84"), pch = 19, cex = 0.3, add = TRUE, col = "red")

```

Otherwise you could use geosphere::destPoint to project out to a certain distance from a centre and get lon/lats for the circle, then transform those to Mercator.

mdsumner
  • 8,196
  • 1
  • 27
  • 30
  • Thank you for your answer @mdsummer. However, I don't understand why you have to reproject. The main idea is to write a parametric equation of a circle with lat lon, and the output is also lat lon. when I simply draw these coordinates on a map (which is done with a mercator projection), we then can see the form of the circle on the map. For me, the projection is done with leaflet. – John Smith Jun 01 '16 at 10:11
  • I missed this comment. It depends what you mean by a "circle". I'm thinking of a radius in real-world distance from a centre point. One way or another you need to find what "equi-distant" from that point is in a longitude/latitude or Mercator or whatever coordinate system. If you mean equi-distant in planar coords in long-lat, that's different but I suspect you don't mean that. – mdsumner May 22 '17 at 04:52