grids in the right place

gridstuff.R

# make a boundary

library(mapdata)
library(maptools)
library(rgeos)
library(ggplot2)
library(sp)

# build the coastline
coastline <- map("world", c("USA", "Canada"), plot=FALSE)
coastline <- map2SpatialPolygons(coastline, ID=coastline$name,
                                 proj4string=CRS("+proj=longlat +datum=WGS84"))

# box around the area (change these to get a smaller area)
outer <- cbind(c(-69.5, -69.5, -71, -71, -69.5),
               c(41.5, 42.5, 42.5, 41.5, 41.5))
outer <- Polygons(list(region=Polygon(coords=outer)), ID="outer")
outer <- SpatialPolygons(list(outer),
                         proj4string=CRS("+proj=longlat +datum=WGS84"))

# get the intersection (well, difference)
gg <- gIntersection(outer, coastline)
bnd <- gSymdifference(outer, gg)

# check what it looks like
plot(bnd)

# make some points
# use this to change # points
res <- 200
dat <- expand.grid(x=seq(min(outer[, 1]), max(outer[, 1]), len=res),
                   y=seq(min(outer[, 2]), max(outer[, 2]), len=res))

# get the points inside that boundary (faff)
x <- dat$x
y <- dat$y
library(mgcv)

fbnd <- fortify(bnd)
fbnd <- fbnd[fbnd$piece==1, 1:2]
names(fbnd) <- c("x", "y")
ind <- inSide(fbnd, x, y)

# check what we got
plot(dat[ind, ])
plot(bnd, add=TRUE)