5

Given a polygon in a PostGIS database, what is the best way to create a graph of its distribution over a range of latitudes? I'm picturing a graph of a curve, where the proportion of the total area under the curve that lies between two points on the x-axis (latitude) represents the proportion of the polygon's area between those latitudes. The tools I have at my disposal are PostGIS and R. It seems like this should be a fairly routine exercise for which there are well established approaches that I simply haven't been able to find.

Gregory
  • 1,661
  • 4
  • 20
  • 30
  • Create a number of horizontal (or vertical?) slice rectangles that span your polygon. Give each one a latitude attribute. Then do an ST_overlay or ST_union or something, so that your polygon is sliced. Then compute the area of each slice. Would that do it? I've put this in a comment due to vagueness... – Spacedman May 01 '12 at 21:49
  • A simple modification of the line-sweep algorithm at http://gis.stackexchange.com/a/33449/664 will do it: you only have to compute the length of the intersection of the sweep line with the polygon at each event, provided you use an cylindrical equal-area projection. The result will be the density function; integrating that will give the requested cumulative function. – whuber Jun 05 '15 at 16:33

1 Answers1

4

OKay, here in an answer, some R code - uses rgeos, sp:

For p = a single SpatialPolygon or single row of a SpatialPolygonsDataFrame:

slice <- function(p,n=20){
  bb = bbox(p)
  ys = seq(bb[2,1],bb[2,2],len=n)
  ll = list()
  for(s in 1:(n-1)){
    ll[[s]] = Polygons(list(
        Polygon(cbind(
                      c(bb[1,1],bb[1,2],bb[1,2],bb[1,1],bb[1,1]),
                      c(ys[s],ys[s],ys[s+1],ys[s+1],ys[s])))),s)
  }
  SpatialPolygons(ll)
}

where n=number of slices, returns n rectangular slice polygons. Now overlay:

lls = slice(p,40)
ii = gIntersection(p,lls,byid=TRUE)
plot(ii)

should show you the slices. Now get the areas:

a = sapply(slot(ii,"polygons"),slot,"area")

I'm just not sure how correct that will be if the polygons have holes in them. Try some tests first.

Spacedman
  • 63,755
  • 5
  • 81
  • 115
  • Thanks, that looks like it might work. I just installed rgeos and I'll post back when I have a chance to try it out, unfortunately it might be a few days or more. – Gregory May 02 '12 at 00:47