pi

library(sp)

## Warning: package 'sp' was built under R version 3.2.5

library(raster)

## Warning: package 'raster' was built under R version 3.2.5

library(rgdal)

## Warning: package 'rgdal' was built under R version 3.2.5

## rgdal: version: 1.2-5, (SVN revision 648)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.1.2, released 2016/10/24
##  Path to GDAL shared files: 
##  Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
##  Path to PROJ.4 shared files: (autodetected)
## WARNING: no proj_defs.dat in PROJ.4 shared files
##  Linking to sp version: 1.2-3

library(rgeos)

## Warning: package 'rgeos' was built under R version 3.2.5

## rgeos version: 0.3-22, (SVN revision 544)
##  GEOS runtime version: 3.4.2-CAPI-1.8.2 r3921 
##  Linking to sp version: 1.2-4 
##  Polygon checking: TRUE

library(plyr)

## Warning: package 'plyr' was built under R version 3.2.5

library(dplyr)

## Warning: package 'dplyr' was built under R version 3.2.5

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:plyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize

## The following objects are masked from 'package:rgeos':
## 
##     intersect, setdiff, union

## The following objects are masked from 'package:raster':
## 
##     intersect, select, union

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(ggplot2)

## Warning: package 'ggplot2' was built under R version 3.2.5

library(scales)

## Warning: package 'scales' was built under R version 3.2.5

set.seed(1)

square <- readWKT('POLYGON((-1 -1,-1 1,1 1,1 -1, -1 -1))')
circle <- readWKT('POINT(0 0)') %>% 
  gBuffer(width = 1, quadsegs = 25)
plot(square, axes = F, border = '#FA6900', lwd = 2)
plot(circle, add = T, col = '#69D2E7', border = 'transparent')

n_pts <- 1000
pts <- spsample(square, n_pts, type = 'random')
plot(square, axes = F, border = '#FA6900', lwd = 2)
plot(circle, add = T, col = '#69D2E7', border = 'transparent')
plot(pts, add = T, pch = 21, cex = 0.25, col = '#333333')

over(pts, circle) %>% 
  {!is.na(.)} %>% 
  sum

## [1] 768

(pts_within <- pts[circle, ])

## class       : SpatialPoints 
## features    : 768 
## extent      : -0.9738448, 0.9855762, -0.9905714, 0.981201  (xmin, xmax, ymin, ymax)
## coord. ref. : NA

length(pts_within)

## [1] 768

n_within <- gIntersects(pts, circle, byid = T) %>% 
  sum
n_within

## [1] 768

(pi_est <- 4 * n_within / n_pts)

## [1] 3.072

round(abs(100 * (pi_est / pi - 1)), 2)

## [1] 2.22

estimate_pi <- function(n_pts) {
  t <- system.time({
    square <- readWKT('POLYGON((-1 -1,-1 1,1 1,1 -1, -1 -1))')
    circle <- readWKT('POINT(0 0)') %>% 
      gBuffer(width = 1, quadsegs = 25)
    n_within <- spsample(square, n_pts, type = 'random') %>% 
      gIntersects(circle, byid = T) %>% 
      sum
    pi_estimate <- (4 * n_within / n_pts)
  })
  return(data.frame(pi_estimate, t = t[['user.self']]))
}
estimate_pi(1e4)

##   pi_estimate     t
## 1      3.1408 0.322