Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I am looking automate the creation of a polar grid shapefile (polygon). I am currently creating this in a many-step editing process, starting from a concentric buffer.

Here is how I think I should be approaching the issue. Given a polar grid of this shape:

enter image description here

Where V1 represents the origin, and the length of the line segment (v1,v4) represents the buffer distance. V2 and V3 represent points along the concentric ring that intersects the polar lines. The amount of polar lines can be changed.

Any thoughts on how to do this?

share|improve this question
    
What software are you using (or can use)? –  whuber Aug 10 '12 at 23:02
    
Going to need to run my analysis in ArcMap so the final product will be a shapefile/featureclass. I am fairly proficient in Python and have experience with R. –  Michael Markieta Aug 11 '12 at 2:29
    
Also add Mathematica, but I am not sure if that will aid at all (there are polar grid functions). –  Michael Markieta Aug 11 '12 at 2:35
add comment

1 Answer

up vote 4 down vote accepted

It's not really any different than creating a Cartesian grid; the only complication involves representing the circular arcs that bound each cell.

As in the Cartesian case, the input should include an origin, a set of radii to use, and a set of angles to use. For greater generality, we needn't require that radii or angles be equally spaced.

Here's an example in R showing a grid with varying radial widths (but not starting at the origin!) and constant angular widths (but not necessarily covering a full circle). It represents the cells as polygons, because this is useful for typical grid applications, such as collecting points within cells (a point-in-polygon operation).

require(shapefiles)
shape.polygon <- 5 # Use 3 for polyline
#
# Grid specification.
#
origin <- c(100,200)
r <- c(1,3,6,10)                                       # In same units as the origin coordinates
theta <- seq(from=-80, to=150, length.out=14) * pi/180 # Must be in radians
#
# Create and write the grid.
#
g <- radial.grid(r, theta, origin, 6 * pi/180)
shapefile <- convert.to.shapefile(g$shapes, g$attrs, "Id", shape.polygon)
write.shapefile(shapefile, "f:/temp/grid", arcgis=T)

Example

Obviously the hard work is done by radial.grid. The boundary of a cell is traversed by following its outer arc, jumping to the end of the inner arc, following that in reverse order, and then returning to the beginning of the first arc. Polygonal approximations to these arcs create vertices approximately every e radians (controlled by the caller). The work of creating these arcs is halved by computing one of them and then rescaling and reversing it to create the other. (The work could be halved again by not recomputing arcs shared by cells, but I didn't bother to make this optimization.)

Within radial.grid, then, there are two helper functions: arc creates a unit-radius arc from one angle to another. (Angles are measured with the mathematical convention: in radians counterclockwise from due east. To measure [still in radians] clockwise from north, just switch sin and cos.) cell uses arc as described to create one polygonal cell. The iteration over all cells is performed by lapply to create a list of polygons. The shapefiles package needs this list to be flattened into a matrix. Finally, radial.grid creates some useful attributes: a sector index, a radial index, and a unique integral identifier.

radial.grid <- function(r, theta=c(0, 2*pi), origin=c(0,0), e=(pi/180)) {
  if (length(r) < 2) stop("Provide at least two radii.")
  if (length(theta) < 2) stop("Provide at least two angles.")
  # origin: (x,y) coordinates of center of grid
  # r:      vector of radii
  # theta:  vector of angles
  # e:      angular increment to approximate circular arcs
  #
  # Returns a list of dataframes convertible to a shapefile, named 
  # `shapes` and `attrs`.
  #-----------------------------------------------------------------------------#
  #
  # Create a sequence of cartesian coordinates, as a 2 by m matrix, tracking
  # in the positive direction along a circular arc from (1,theta1) to (1,theta2).
  #
  arc <- function(theta1, theta2) {
    dtheta <- (theta2 - theta1) %% (2 * pi)  # Assure a positive orientation
    theta <- seq(from=theta1, to=theta1+dtheta, length.out=ceiling(dtheta/e)+1)
    sapply(theta, function(a) c(cos(a), sin(a)))
  }
  #
  # Create a polygon with identifier `id` for a radial cell bounded by angles 
  # `theta1` and `theta2` and radii `r1` < `r2`.
  #
  cell <- function(id, origin=c(0,0), r1, r2, theta1, theta2) {
    a <- arc(theta1, theta2)
    if (r1==0) b <- matrix(c(0,0), nrow=2)
    else b <- r1 * a[, dim(a)[2]:1]
    a <- r2 * a
    xy <- cbind(a, b, a[,1])
    rbind(rep(id, dim(xy)[2]), xy + origin)
  }
  #
  # Preprocess the input: make sure radii and angles are sorted.
  #
  rho <- sort(r)
  th <- sort(theta)
  #
  # Create the cells delimited by `rho` and `th`.
  #
  n.annuli <- (length(rho)-1) 
  n.wedges <- (length(th)-1)
  ixy <- lapply(1:(n.wedges * n.annuli)-1, function(ij) {
                  i <- ij %% n.wedges + 1; j <- floor(ij/n.wedges)+1;
                  cell(ij+1, origin, rho[j], rho[j+1], th[i], th[i+1])
                })           # List of identified polygons
  #
  # Prepare the cells for writing in shapefile format.
  #
  ixy <- do.call(cbind, ixy) # Flatten into a matrix
  df.shapes <- data.frame(Id=ixy[1,], X=ixy[2,], Y=ixy[3,])
  df.table <- expand.grid(1:n.wedges, 1:n.annuli)
  colnames(df.table) <- c("Sector", "Radius")
  df.table$Id <- df.table$Sector + n.wedges*(df.table$Radius-1)
  list(shapes=df.shapes, attrs=df.table)
}
share|improve this answer
    
Fascinating stuff. Thank you @whuber –  Michael Markieta Aug 13 '12 at 13:08
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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