Clip a big number of parcels to the University District. All data downloaded from wagda.lib.washington.edu (must have netid to download)

big-clip.R

# install.packages(c("rgdal", "rgeos", "ggmap", "maptools"))
library(rgdal)
library(rgeos)
library(ggmap)
library(maptools)

# Read the shapefile with OGR. Punctuation matters: No trailing slash
# on the dsn, and no .shp on the layer
udist <- readOGR(dsn="/Users/peter/Downloads/clip",
                 layer="unidstuc")

# Common projection: WGS84
udist <- spTransform(udist, CRS("+init=EPSG:4326"))

# This one takes an obnoxiously long time. Still < 3 min on my machine
rparcels <- readOGR(dsn="/Users/peter/Downloads/clip",
                    layer="parcel")

# Common projection: WGS84, first setting the (possibly) correct
# string from http://spatialreference.org/ref/epsg/ Also obnoxiously
# slow. ~ 3 min.

rparcels@proj4string <- CRS("+proj=lcc +lat_1=48.73333333333333 +lat_2=47.5 +lat_0=47 +lon_0=-120.8333333333333 +x_0=500000.0001016001 +y_0=0 +ellps=GRS80 +datum=NAD83 +to_meter=0.3048006096012192 +no_defs")

rparcels <- spTransform(rparcels, CRS("+init=EPSG:4326"))

# Return matrix that determines whether features in rparcels layer
# intersects udist layer. Note, if you had more than one feature in
# udist, you'll need to flatten the resulting matrix horizontally
# before converting to vector.
i <- as.vector(gIntersects(udist, rparcels, byid = TRUE))

# Technically not a clip yet. Parcels are kept whole if they touch the
# UD.
uparcels <- rparcels[i,]

# The clip runs much more quickly on the small set, though you might
# not actually want to be clipping...

# uparcels <- gIntersection(udist, uparcels, byid = TRUE)
prettymap <- ggmap(get_map(location=uparcels@bbox, maptype='toner'))

# ggplot and ggmap like data to be in dataframe format, not SPDF object
uparcels.df <- fortify(uparcels)

# Join the data back in. What a pain!
uparcels@data$id <- row.names(uparcels) # common id
uparcels.df <- merge(x=uparcels.df, y=uparcels@data, by="id")

prettymap + geom_polygon(data=uparcels.df,
                         mapping=aes(x=long, y=lat, group=group, fill=PU_CAT_DES),
                         alpha=0.5)