Polar grids in ArcGIS

Question

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:

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?

Answer 1

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)

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)
}