Dici che il fiume trova la via al mare e come il fiume giungerai a me (Miss Sarajevo, U2)

One way to calculate approximately the area of some place is to circumscribe it into a polygon of which you know its area. After that, generate coordinates inside the polygon and count how many of them fall into the place. The percentage of coordinates inside the place by the area of the polygon is an approximation of the desired area.

I applied this technique to calculate the area of the Vatican City. I generated a squared grid of coordinates around the Capella Sistina (located inside the Vatican City). To calculate the area I easily obtain the convex hull polygon of the coordinates using `chull`

function of `grDevices`

package. Then, I calculate the area of the polygon using `areaPolygon`

function of `geosphere`

package.

To obtain how many coordinates of the grid *fall* inside the Vatican City, I use `revgeocode`

function of `ggmap`

package (I love this function). For me, one coordinate is inside the Vatican City if its related address contains the words “Vatican City”.

What happens generating a grid of 20×20 coordinates? I obtain that the area of the Vatican City is about 0.32Km^{2 } but according to Wikipedia, the area is 0.44Km^{2}: this method underestimates the area around a 27%. But why? Look at this:

This plot shows which addresses of the grid fall inside the Vatican City (ones) and which of them do not fall inside (zeros). As you can see, there is a big zone in the South, and a smaller one in the North of the city where reverse geocode do not return “Vatican City” addresses.

Maybe Pope Francis should phone Larry Page and Sergey Brin to claim this 27% of his wonderful country.

I was willing to do this experiment since I wrote this post. This is the code:

require(geosphere) require(ggmap) require(plotGoogleMaps) require(grDevices) setwd("YOUR-WORKING-DIRECTORY-HERE") #Coordinates of Capella Sistina capella=geocode("capella sistina, Vatican City, Roma") #20x20 grid of coordinates around the Capella g=expand.grid(lon = seq(capella$lon-0.010, capella$lon+0.010, length.out=20), lat = seq(capella$lat-0.005, capella$lat+0.005, length.out=20)) #Hull Polygon containing coordinates p=g[c(chull(g),chull(g)[1]),] #Address of each coordinate of grid a=apply(g, 1, revgeocode) #Estimated area of the vatican city length(grep("Vatican City", a))/length(a)*areaPolygon(p)/1000/1000 s=cbind(g, a) s$InOut=apply(s, 1, function(x) grepl('Vatican City', x[3]))+0 coordinates(s)=~lon+lat proj4string(s)=CRS('+proj=longlat +datum=WGS84') ic=iconlabels(s$InOut, height=12) plotGoogleMaps(s, iconMarker=ic, mapTypeId="ROADMAP", legend=FALSE)