0

Similarly to my previous question here, say I am plotting countries on a world map using maptools, if I were to plot a country, is there a way of finding the central point of this country and plotting a radial distance from this point? I am using the shapefile wrld_simpl that comes with maptools, so if I plot Germany:

 plot(wrld_simpl[wrld_simpl$NAME=='Germany',], col='red', add=T)

I would want to find the centre of Germany and plot a circle from this point that shows any area that falls within say a 100km distance. Again I want to be able to do this for lots of different countries so I'd ideally want a general solution, not one specific to just Germany.

Community
  • 1
  • 1
userk
  • 901
  • 5
  • 11
  • 19

1 Answers1

1

This is a partial answer. The code snippet below calculates the centroid of a polygon, so if you can pull the polygon vertex data for your country of interest, this will give you the "center," after which it's trivial to draw a circle. (polyx and polyy are vectors of x- and y- coordinates

require(pracma)
pchit <- polyarea(polyx,polyy)
centx <- centy <- 0
    for (kk in 1:(length(polyx)-1) ) {
        centx <- centx + (polyx[kk]+polyx[kk+1]) * (polyx[kk]*polyy[kk+1]-polyx[kk+1]*polyy[kk])
        centy <- centy + (polyy[kk]+polyy[kk+1]) * (polyx[kk]*polyy[kk+1]-polyx[kk+1]*polyy[kk])
    }
    centx <- -1/pchit/6 * centx
    centy <- -1/pchit/6 * centy
Carl Witthoft
  • 20,573
  • 9
  • 43
  • 73
  • Thanks for the tip. I will add a vectorized form of this function as polycenter() to the package. The `-1`seems to be wrong. Is there any reference for this formula? – Hans W. Aug 22 '13 at 09:03
  • I changed `polyarea` that now returns the true (oriented) area, not its absolute value. This implies the `-1` is not necessary any more. – Hans W. Aug 22 '13 at 12:31
  • @HansWerner Apologies - you meant the minus sign in the last two lines, right? Yes, those go away depending on the ordering of the vertices. – Carl Witthoft Aug 22 '13 at 13:10