# Create spatial polygon grid from spatial points in R

I am trying to recreate a map in R. I need 2 by 5 degree grid cells across the contiguous United States, but plotted in Albers Equal Area. The result is a grid with curved lines and slightly variable grid size, so I assume I cannot use a raster. I am able to reproduce the points, but not the grid. How do I make grid (polygon) lines for these points?

``````library(raster)

e <- as(raster::extent(-125, -65, 25, 49), "SpatialPolygons")
grd_lrg <- as.data.frame(makegrid(e, "regular", cellsize = c(5,2)))
names(grd_lrg)       <- c("X", "Y")
coordinates(grd_lrg) <- c("X", "Y")

# add projection information to the empty grid
proj4string(grd_lrg) <- CRS("+init=epsg:4326")

# transform to equal area Albers USGS
grd_lrg <- spTransform(grd_lrg, CRS("+init=epsg:5070"))

# plot it
plot(grd_lrg)
`````` If I create the grid in Arc Map using Create Fishnet (Data Management) with the WGS84 coordinate system and import it into R, I get a nice looking grid: But converting it to Albers Equal Area with spTransform creates major distortions of the lines: So I am thinking I need to convert the points to Albers first and then add the lines?

• That distortion appears to be because each line across the fishnet is defined only by its endpoints, so the transformed line is going to be a straight line between the transformed endpoints, rather than curved. Is there an option in the fishnet generator to make polygons rather than lines? Then each square would have four points that would transform better. Mar 2, 2018 at 12:51
• Yes, when I selected the polygon option (polyline is default in Arc Map create fishnet) and then transformed to Albers the result is good. Mar 2, 2018 at 17:41

here is a solution using sf and sf::st_make_grid:

``````library(raster)
library(sf)

e <- as(raster::extent(-125, -65, 25, 49), "SpatialPolygons") %>%
st_as_sf()

grd_lrg <- st_make_grid(e, cellsize = c(5, 2)) %>%
st_set_crs(4326) %>%
st_transform(5070)

plot(grd_lrg)
`````` • perhaps add 'library (dplyr)' to this answer? Aug 18, 2020 at 10:56

You can make a long/lat raster and transform

``````library(raster)
r <- raster(ext = extent(-125, -65, 25, 49), res=c(5,2))
values(r) <- 1:ncell(r)
rA <- projectRaster(r, crs="+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m")
``````

But what you want seems to be

``````library(rgdal)
p <- rasterToPolygons(r)
# or do:  p <- as(r, "SpatialPolygonsDataFrame")
pA <- spTransform(p, "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=23 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m")

plot(rA, axes=FALSE, box=FALSE) • You can interpolate to points, e.g. `centroids(pA)` . Or to a raster and then use `extract(r, pA, fun=mean)` . If that is not clear, you may want to set up another question with some example data. Mar 2, 2018 at 1:14