I am attempting to plot a colored grid of attributes over a base map. Similar questions have been asked here, including one recent question I asked myself, but earlier solutions have usually involved ggplot2. I am hoping to find a solution ideally without using ggplot2 if only because my final actual map will be quite complex and I have made substantial progress so far using plot.
Below is code that plots Colorado and creates a grid of fake attribute data. I can convert the grid into spatial points using SpatialPoints and plot those points. However, I suspect I need to convert to spatial polygons, perhaps using SpatialPolygons. When I attempt to use SpatialPolygons I obtain the error shown below.
Once I get the grid cell layer added to the map I would like to crop the grid to only be displayed within Colorado (or, for example, within Colorado and Wyoming if I added another state to the map and used a larger attribute grid). However, that perhaps is a follow-up question.
Here is an earlier question with an answer using plot, but requesting additional information from the poster: R Plot Filled Longitude-Latitude Grid Cells on Map
Here is a similar earlier question from myself that only has a ggplot2 solution: color grid cells in the United States and Canada
My approach below is based off of a solution to another of my earlier questions: R: creating a map of selected Canadian provinces and U.S. states
Thank you for any help.
library(rgdal)
library(maptools)
library(ggplot2)
library(plyr)
library(RColorBrewer)
library(classInt)
library(raster)
# NW(long,lat) SE(long,lat)
mapExtent <- rbind(c(-115, 43), c( -97, 35))
# assign projection
newProj <- CRS("+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0")
## Project map extent
mapExtentPr <- spTransform(SpatialPoints(mapExtent,
proj4string=CRS("+proj=longlat")), newProj)
us1 <- getData('GADM', country="USA", level=1)
colorado <- us1[(us1$NAME_1 %in% c('Colorado')),]
## Project Colorada layer
coloradoPr <- spTransform( colorado, newProj)
# Create grid cells containing fake attribute data
# using approach found here:
# http://www.numbertheory.nl/2011/11/08/drawing-polar-centered-spatial-maps-using-ggplot2/
set.seed(1234)
xlim = c(-113, -99)
ylim = c( 35, 45)
dat.grid = expand.grid(x = xlim[1]:xlim[2], y = ylim[1]:ylim[2])
dat.grid$z = runif(nrow(dat.grid))
## Project grid layer
dat.gridPr <- spTransform(SpatialPoints( dat.grid, proj4string=CRS("+proj=longlat")), newProj)
dat.gridPr2 <- spTransform(SpatialPolygons(dat.grid, proj4string=CRS("+proj=longlat")), newProj)
#Error in spTransform(SpatialPolygons(dat.grid, proj4string = CRS("+proj=longlat")), :
# error in evaluating the argument 'x' in selecting a method for function 'spTransform': Error in SpatialPolygons(dat.grid, proj4string = CRS("+proj=longlat")) :
# cannot get a slot ("area") from an object of type "integer"
## Plot each projected layer, beginning with the projected extent
plot(mapExtentPr, pch=NA)
plot(coloradoPr , border="white", col="lightgrey", add=TRUE)
plot(dat.gridPr , border="white", col="lightgrey", add=TRUE)
plot(dat.gridPr , fill=dat.grid$z, add=TRUE)
plot(dat.gridPr , fill=dat.gridPr$z, add=TRUE)
EDIT
This post describes how to add grid cells one at a time, which I could do if need be:
