31

How can one use R to

  1. split a shapefile in 200 meter squares/sub-polygons,
  2. plot this grid (incl. ID numbers for each square) over the original map below, and
  3. evaluate in which square specific geographic coordinates are located.

I am a beginner in GIS and this is perhaps a basic question, but I haven't found a tutorial on how to do this in R.

What I have done so far is loading a shapefile of NYC and plotting some exemplary geographic coordinates.

I appreciate any general advice, but I am looking for an example (R code) how to this with the data below.

# Load packages 
library(maptools)

# Download shapefile for NYC
# OLD URL (no longer working)
# shpurl <- "http://www.nyc.gov/html/dcp/download/bytes/nybb_13a.zip"
shpurl <- "https://www1.nyc.gov/assets/planning/download/zip/data-maps/open-data/nybb_13a.zip"

tmp    <- tempfile(fileext=".zip")
download.file(shpurl, destfile=tmp)
files <- unzip(tmp, exdir=getwd())

# Load & plot shapefile
shp <- readShapePoly(files[grep(".shp$", files)])
plot(shp)

# Define coordinates 
points_of_interest <- data.frame(y=c(919500, 959500, 1019500, 1049500, 1029500, 989500), 
                 x =c(130600, 150600, 180600, 198000, 248000, 218000),
                 id  =c("A"), stringsAsFactors=F)

# Plot coordinates
points(points_of_interest$y, points_of_interest$x, pch=19, col="red")

enter image description here

36
+100

Here is an example using a SpatialGrid object:

### read shapefile
library("rgdal")
shp <- readOGR("nybb_13a", "nybb")

proj4string(shp)  # units us-ft
# [1] "+proj=lcc +lat_1=40.66666666666666 +lat_2=41.03333333333333 
# +lat_0=40.16666666666666 +lon_0=-74 +x_0=300000 +y_0=0 +datum=NAD83
# +units=us-ft +no_defs +ellps=GRS80 +towgs84=0,0,0"

### define coordinates and convert to SpatialPointsDataFrame
poi <- data.frame(x=c(919500, 959500, 1019500, 1049500, 1029500, 989500),
                  y=c(130600, 150600, 180600, 198000, 248000, 218000),
                  id="A", stringsAsFactors=F)
coordinates(poi) <- ~ x + y
proj4string(poi) <- proj4string(shp)

### define SpatialGrid object
bb <- bbox(shp)
cs <- c(3.28084, 3.28084)*6000  # cell size 6km x 6km (for illustration)
                                # 1 ft = 3.28084 m
cc <- bb[, 1] + (cs/2)  # cell offset
cd <- ceiling(diff(t(bb))/cs)  # number of cells per direction
grd <- GridTopology(cellcentre.offset=cc, cellsize=cs, cells.dim=cd)
grd
# cellcentre.offset 923018 129964
# cellsize           19685  19685
# cells.dim              8      8

sp_grd <- SpatialGridDataFrame(grd,
                               data=data.frame(id=1:prod(cd)),
                               proj4string=CRS(proj4string(shp)))
summary(sp_grd)
# Object of class SpatialGridDataFrame
# Coordinates:
#      min     max
# x 913175 1070655
# y 120122  277602
# Is projected: TRUE
# ...

Now you can use the implemented over-method to obtain the cell IDs:

over(poi, sp_grd)
#   id
# 1 57
# 2 51
# 3 38
# 4 39
# 5 14
# 6 28

To plot the shapefile and the grid with the cell IDs:

library("lattice")
spplot(sp_grd, "id",
       panel = function(...) {
         panel.gridplot(..., border="black")
         sp.polygons(shp)
         sp.points(poi, cex=1.5)
         panel.text(...)
       })

spplot1

or without colour/colour key:

library("lattice")
spplot(sp_grd, "id", colorkey=FALSE,
       panel = function(...) {
         panel.gridplot(..., border="black", col.regions="white")
         sp.polygons(shp)
         sp.points(poi, cex=1.5)
         panel.text(..., col="red")
       })

spplot2

  • This looks like an answer to me, but in case you are looking for something different. Try the r tag in stackoverflow stackoverflow.com/search?q=R+tag – Brad Nesom Mar 14 '14 at 2:26
  • @rcs this code looks just like what I am trying to do but my shapefile is in a different projection: proj4string (DK_reg1) [1] "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" does anyone have any suggestions on how to break this shapefiles of this projection to 1000 equal sized grid cells? and then randomly select 100 of them and highlight them? – I Del Toro Sep 13 '14 at 11:34
8

The New York dataset provided in the question is no longer available for download. I use the nc dataset from sf package to demonstrate a solution using sf package:

library(sf)
library(ggplot2)

# read nc polygon data and transform to UTM 
nc <- st_read(system.file('shape/nc.shp', package = 'sf')) %>%
  st_transform(32617)

# random sample of 5 points
pts <- st_sample(nc, size = 5) %>% st_sf

# create 50km grid - here you can substitute 200 for 50000
grid_50 <- st_make_grid(nc, cellsize = c(50000, 50000)) %>% 
  st_sf(grid_id = 1:length(.))

# create labels for each grid_id
grid_lab <- st_centroid(grid_50) %>% cbind(st_coordinates(.))

# view the sampled points, polygons and grid
ggplot() +
  geom_sf(data = nc, fill = 'white', lwd = 0.05) +
  geom_sf(data = pts, color = 'red', size = 1.7) + 
  geom_sf(data = grid_50, fill = 'transparent', lwd = 0.3) +
  geom_text(data = grid_lab, aes(x = X, y = Y, label = grid_id), size = 2) +
  coord_sf(datum = NA)  +
  labs(x = "") +
  labs(y = "")

# which grid square is each point in?
pts %>% st_join(grid_50, join = st_intersects) %>% as.data.frame

#>   grid_id                 geometry
#> 1      55 POINT (359040.7 3925435)
#> 2      96   POINT (717024 4007464)
#> 3      91 POINT (478906.6 4037308)
#> 4      40 POINT (449671.6 3901418)
#> 5      30 POINT (808971.4 3830231)

enter image description here

  • Thanks. I updated the link in my question to relfect the changes on their webpage. Now it should work again. – majom Mar 5 '18 at 15:18
  • I really need to start using the sf package. This is awesome! – spacedSparking Mar 23 '18 at 20:41
  • Is there an easy way to only plot the grid cells that intersect with the state polygon? – spacedSparking Mar 23 '18 at 22:09
  • st_intersection(grid_50, nc) should do it – sebdalgarno Mar 24 '18 at 23:18
2

Have you looked at the R raster package? It has tools to convert to/from vector GIS objects so you should be able to a) create a raster (grid) with 200x200m cells and b) convert it to a set of polygons with a logical id of some kind. From there I would look at the sp package to help with intersecting the points and the polygon grid. This http://cran.r-project.org/web/packages/sp/vignettes/over.pdf page might be a good start. Wandering through the sp package docs you might be able to start with the SpatialGrid-class and just skip the raster part entirely.

  • 1
    Thanks for your comment. I already looked at those packages, but failed to implement this. This was why I was asking for an example. – majom Mar 10 '14 at 9:05
0

The "GIS universe" is complex and have many standards that your data must be compliant. All "GIS tools" interoperates by GIS-standards. All "serious GIS data" today (2014) are stored in a database.

The best way to "use R" in a GIS context, with other FOSS tools, is embedded into SQL. The best tools are PostgreSQL 9.X (see PL/R) and PostGIS.


You answer:

  • To import/export shape files: use shp2pgsql and pgsql2shp.
  • To "split a shape file in 200 meter squares/sub-polygons": see ST_SnapToGrid(), ST_AsRaster(), etc. We need understand better your needs to express into a "recipe".
  • you say that need "geographic coordinates are located" .. perhaps ST_Centroid() of the squares (?)... You can express "more mathematically" so I understand.

... Perhaps you not need any raster convertion, only a matrix of regurlar-sampled points.


A primitive way is use R without PL/R, in a your usual external compiler: only convert your polygons and export as shape or as WKT (see ST_AsText), then convert data with awk or another filter to the R format.

  • 1
    Thanks for your help. However, I would strongly prefer a solution which relies completely on R and existing packages. When I am able to split the shape file in 200m*200m subpolygons I can check with point.in.polygon which coordinates are in which polygons. My problem is to split the original shapefile in those sub-polygons. – majom Mar 9 '14 at 21:38
  • It is not a "SIG problem", try to ask at Stackoverflow directly. – Peter Krauss Mar 10 '14 at 11:18

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.