2

I have a non-overlapping polygon-based shapefile (.shp) with a large spatial extent and many dozens of associated attributes. The shapefile is projected in UTMs. I would like to convert the polygons to points spaced out in a 30-m resolution grid, in which each point would retain the attributes of the polygon it is located within.

Output would simply be a table of the points: 

X, Y, attribute1, attribute2, attribute 3,etc...

I would ideally like to do this operation in R, or (less ideally) some other free program I can run on a Mac.

Luke
  • 4,769
  • 5
  • 28
  • 36
  • 1
    just an idea, how about creating the 30m resolution grid (as points) to cover the polygon and do a `point.in.polygon` to pull out those that lie on the line of the polygon only. If you have any attributes, you can then attach them to a separate df? – RJ- Feb 08 '13 at 16:54
  • have you tried to attempt anything yourself ? – CHP Feb 08 '13 at 16:58

1 Answers1

5

NOTE: I'm throwing this up in part to learn whether there's a more elegant way to do any of this. So, please, spatial types, pitch in with any suggestions for improvement.

(In particular, Step 2, which sets up a "SpatialPoints" grid with the points to which values will be extracted, always seems painfully low-level to me.)

This uses over() to extract attributes from a "SpatialPolygonDataFrame" at the coordinates contained in a "SpatialPoints" object constructed for just that purpose.

library(rgdal)

## (1) Read in an example shapefile
dsn <- system.file("vectors", package = "rgdal")[1]
scot_BNG <- readOGR(dsn=dsn, layer="scot_BNG")
scot_BNG <- scot_BNG[1:5,]  # Let's just use part of it

## (2) Set up a SpatialPoints object with the grid of points 
##     for which you want to extract values
res <- 10000            ## Distance between grid points (30 in OP's question) 
BB <- bbox(scot_BNG)
BB <- res*round(BB/res) ## Pretty up the bounding box
GT <- GridTopology(cellcentre.offset = BB[,1], 
                   cellsize = c(res, res),
                   cells.dim = (c(diff(BB[1,]), diff(BB[2,]))/res) + 1)
SP <- SpatialPoints(GT, proj4string = CRS(proj4string(scot_BNG)))

## (3) Extract the values
vals <- over(SP, scot_BNG)
res <- cbind(coordinates(SP), vals)

## Finally, have a look at a few of the points.
x <- res[!is.na(res$SP_ID),]
rbind(head(x,3), tail(x,3))[1:10]
#          x      y SP_ID       NAME ID_x COUNT   SMR  LONG  LAT    PY
# 4   230000 970000     0 Sutherland   12     5 279.3 58.06 4.64 37521
# 5   240000 970000     0 Sutherland   12     5 279.3 58.06 4.64 37521
# 25  220000 960000     0 Sutherland   12     5 279.3 58.06 4.64 37521
# 425 260000 780000     4   Bedenoch   17     2 186.9 57.06 4.09 27075
# 426 270000 780000     4   Bedenoch   17     2 186.9 57.06 4.09 27075
# 427 280000 780000     4   Bedenoch   17     2 186.9 57.06 4.09 27075
Josh O'Brien
  • 159,210
  • 26
  • 366
  • 455
  • this should have been solved with the raster package, something like `values(rasterize(myshp, raster(myshp), "attribute_name"))` – Ma Ba Dec 14 '17 at 10:21