9

I want to examine the orientation of each line in a polygon so that I can calculate their solar exposure. Each polygon represents a building, and has an associated height. For the moment, I just want to consider the orientation, and will later consider shading issues.

One approach I thought was to split the polygon into lines, and to calculate the orientation of each line, but the difficulty is that I then have to identify the outer face of that line. Although most of the polygons are simple four-side figures with straight lines, there are a small number where this is not the case (an issue which I just want to consider, but do not have to solve yet).

I'm familiar with python and planned on doing this all from a script.

PolyGeo
  • 65,136
  • 29
  • 109
  • 338
djq
  • 16,297
  • 31
  • 110
  • 182
  • 1
    Do the answers to this question help: http://gis.stackexchange.com/questions/1886/how-to-calculate-orientation-of-line-segments-using-open-source-gis – Sean Nov 22 '10 at 17:52
  • @Sean - thanks, yes that does help in general. I am not familiar with the specific ArcGIS commands that can be used to do this. Also, the issue of knowing what the inside/outside face of the line in a polygon still remains. – djq Nov 22 '10 at 18:39
  • When you say 'orientation', you haven't quite defined it well enough for me --- what's the orientation of a perfect hexagon, for example? – Dan S. Nov 30 '10 at 19:27
  • @Dan S. I just edited my answer to include orientation of each line in the polygon. – djq Dec 01 '10 at 12:24
  • sorry answerers - I didn't manage to assign the bounty properly yesterday. – djq Dec 02 '10 at 13:14

6 Answers6

6

If you just want majority orientation, check out @Mapperz 's answer above.

Otherwise, as you say you could split the polygons into lines using the Polygon To Line tool. This adds a left and right FID field, where the field is -1 if there is no outer polygon - this can cause a bit of mucking about though if your buildings are adjacent or overlap.

From there you could split your lines at every vertex (maybe use Split into COGO lines), and then calculate the angles on each of the lines (potentially by Updating COGO attributes).

Assuming you have your angle field calculated from North the the aspect will be correct where the left_FID is -1, and to get the aspect when the right_FID is -1 just add 180°. Then based on the original FID you can aggregate, get the majority aspect based on length etc.

The Polygon to line tool is scriptable, (as far as I know) COGO tools are not, so you would have to come up with something there yourself.

Hope this helps!

om_henners
  • 15,642
  • 2
  • 46
  • 86
  • 2
    Shameless plug, and a comment: First, if you don't have an Info license, http://code.google.com/p/boundary-generator/ does more or less the same as the Polygon to Line tool. Second -- this gives you a lot more information than just looking at some sort of general orientation of the overall polygon ... for example you can use the resulting table to do much more precise sun exposure calculations, taking into account each individual wall's aspect. – Dan S. Nov 30 '10 at 19:31
5

Calculate Polygon Main Angle (Cartography)

Calculates the dominant angles of input polygon features and assigns the values to a specified field in the feature class

http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//007000000028000000.htm

Mapperz
  • 49,701
  • 9
  • 73
  • 132
5

Finding orientation

Within a script, the polygon will be available as a set of rings--one outer ring and zero or more inner rings--with each ring represented by a cyclically ordered tuple of vectors (v[0], v1, ..., v[m-1], v[m] = v[0]). Each vector gives the coordinates of a vertex (with no two successive vertices coinciding). From this it is straightforward, as others have pointed out, to obtain normal vectors (i.e., vectors perpendicular to the edge directions):

n[i] = t(v[i+1] - v[i]).

The "t" operation rotates a vector 90 degrees counterclockwise:

t( (x, y) ) = (y, -x).

Only the directions of these normal vectors matter, so rescale them to have unit length: a vector (x, y) rescales to (x/s, y/s) where s = Sqrt(x^2 + y^2) (which is the length of its corresponding edge). From now on, let's assume this has been done. Write the components of the resulting unit normal vectors as

n[i] = (u[i], v[i]), i = 0, 1, ..., m-1.

Discriminating outside from inside

As you remark, this leaves a directional ambiguity: should we be using n[i] or -n[i]? Which one points outwards? This question is equivalent to finding the degree of the Gauss map. To compute it, you need to sum the angles by which the normal directions change as you march around a ring. Because the normal vectors have unit length, the cosine of the angle between two successive edges is

Cos(q_i) = n[i] . n[i+1] = u[i]*u[i+1] + v[i]*v[i+1], i = 0, 1, ..., m-1.

(Define n[m] = n[0].)

The sine of the angle between two successive edges is

Sin(q_i) = n[i] . t(n[i+1]) = u[i]*v[i+1] - v[i]*u[i+1].

(Note that these calculations require only sums, differences, and products so far.) Applying the principal inverse tangent function (ATan2) to any such (cosine, sine) pair gives the angle q_i between -180 and 180 degrees. Summing these angles for i = 0, 1, ..., n-1 produces (up to floating point error) the total curvature of the ring, which must be a multiple of 360 degrees; for a closed non-self-intersecting ring it wil be either +360 or -360. In the first case the degree is 1 and in the second case the degree is -1. The normals are all outwardly oriented when the degree of the outer ring is +1 and the degrees of the inner rings are -1. Re-orient them ring by ring, as necessary, according to this rule. That is, if the degree of any ring is the opposite of that which is needed, negate all the normals for that ring. Now you can proceed with your insolation calculations.

whuber
  • 69,783
  • 15
  • 186
  • 281
  • Thank you for a very comprehensive answer. When you say, 'Within a script, the polygon will be available as a set of rings' how is the polygon available? Although familiar with some python scripting I don't know how to interpret the polygon in this way. I am reading through this slowly and trying to understand it, but I have difficulty translating some of the explanation into pseudocode that I can write. – djq Nov 29 '10 at 21:20
  • A few notes on this answer: Typical GIS APIs will always present the outside of a polygon in counterclockwise order, IIRC, which means you don't have to do that fancy discriminating outside from inside -- just use clockwise rotations on the outer rings and counter-clockwise on the holes. A clarification for celenius: the 'set of rings' bit is how most APIs, including ArcGIS's python, allow you to decompose a polygon into its constituent line segments -- a ring is a closed curve of them. Since polygons may support islands and holes, you could have multiple rings... – Dan S. Nov 30 '10 at 19:26
  • @Dan I wish it were the case that GISes maintained such consistent representations. Over the years ESRI software has fluctuated back and forth among negatively and positively oriented polygons and just leaving them unoriented. At this point I'm not sure I would rely even on documented statements concerning their current approach: I would be cautious. It costs little to check the orientation after you have already processed all the edges in a ring. – whuber Nov 30 '10 at 19:31
  • .. well, I did put in that little IIRC disclaimer, thankfully. ;) It's not all that often I do geometry-level stuff on the ESRI stack anymore. – Dan S. Nov 30 '10 at 19:54
  • @Dan S. - I'm still a little unclear about how this can be programmed in python. Are there any guidelines for this? – djq Dec 01 '10 at 12:27
  • @Celenius An example of how the list-of-rings representation is made available, see this answer. Briefly, .PartCount tells you how many rings there are and GetPart() provides each ring as a list of coordinate pairs. See the ESRI help for more. – whuber Jan 07 '12 at 00:47
1

Can this help?

julien
  • 10,186
  • 6
  • 55
  • 93
  • That's an interesting paper you dug up. However, none of the methods described there is relevant to a solar exposure calculation (unless the calculation is intended to be a crude approximation, which does not appear to be the case here). – whuber Nov 27 '10 at 02:09
1

/* Maybe this helps:

Azimuth - pi/2 is the outward facing orientation of the sides of a RHR polygon:

Here is a PostGIS example, you can create the bldg117862 table using the statement at the end. The SRID is EPSG 2271 (PA StatePlane North Feet) and the geometry is a Multipolygon. To visualize in ArcGIS 10, paste query/subqueries into a Query Layer connection to postgis after creating the table bldg117862 .*/

--===START OF QUERY===

/*Outer query provides orientation of outward orthogonals, and creates outward orthogonal lines of equal lengths as those of the sides from the sides' midpoint.

Dominant facing direction(s) will be the sum of the length, grouped by orientation, in descending order */

SELECT line_id as side_id, length, degrees(orthoaz)as orientation,st_makeline(st_setsrid(st_line_interpolate_point(geom,.5),2271), st_setsrid(st_makepoint(st_x(st_line_interpolate_point(geom,.5))+ (length * (sin(orthoaz))), st_y(st_line_interpolate_point(geom,.5))+ (length * (cos(orthoaz)))),2271)) as geom from

--next outer subquery makes lines from sides' point pairs, calculate the azimuth (orthoaz) of outward orthogonal for each segment

( SELECT bldg2009gid, line_id, st_length(st_makeline(startpoint,endpoint))::numeric(10,2) as length, azimuth(startpoint,endpoint),azimuth(startpoint,endpoint) - pi()/2 as orthoaz, st_makeline(startpoint,endpoint) as geom from

/*innermost subquery - use generate_series() to decompose building polygons into sides' startpoint/endpoint point pairs - note1 - force right-hand-rule to ensure common orientation of all polygon sides note2 - example uses multipolygon, for polygon the geometryn() can be removed */

( SELECT generate_series(1, npoints(exteriorring(geometryn(st_forceRHR(geom),1)))-1) as line_id,gid as bldg2009gid, pointn(exteriorring(geometryn(st_forceRHR(geom),1)) ,generate_series(1, npoints(exteriorring(geometryn(st_forceRHR(geom),1)))-1)) as startpoint, pointn(exteriorring(geometryn(st_forceRHR(geom),1)) ,generate_series(2, npoints(exteriorring(geometryn(st_forceRHR(geom),1))))) as endpoint from bldg117862 ) as t1 ) as t2

--===END OF QUERY===

-- the bldg117862 table create/insert statements

SET STANDARD_CONFORMING_STRINGS TO ON; SELECT DropGeometryColumn('','bldg117862','geom'); DROP TABLE "bldg117862"; BEGIN; CREATE TABLE "bldg117862" (gid serial PRIMARY KEY, "motherpin" varchar(14), "taxpin" varchar(14), "status" varchar(15), "area" numeric, "prev_area" numeric, "pct_change" numeric, "picture" varchar(133), "mappage" varchar(6), "sref_gid" int4, "e_address" varchar(19), "a_address" varchar(19), "perim" numeric, "card" int4, "a_addnum" int4, "e_street" varchar(50), "a_street" varchar(50), "e_hsnum" varchar(10)); SELECT AddGeometryColumn('','bldg117862','geom','2271','MULTIPOLYGON',2); INSERT INTO "bldg117862" ("motherpin","taxpin","status","area","prev_area","pct_change","picture","mappage","sref_gid","e_address","a_address","perim","card","a_addnum","e_street","a_street","e_hsnum",geom) VALUES ('33108480001346',NULL,'Changed','2445.59000000','2139.51000000','14.3100000000','http://www.eriecountygov.org/government/assessment/parcelimages.aspx?parcelid=33108480001346" target="_blank">Image Page','33-108','99040','4622 MCCLELLAND AVE','4622 MCCLELLAND AVE','233.37','1','4622','MCCLELLAND AVE',NULL,'4622','0106000020DF080000010000000103000020DF080000010000000B0000008C721D6C98AC34415E2C5BB9D3E32541AE56DE17BEAC34410613E5A0A0E325411AB6C794AEAC3441BA392FE372E32541C89C38429DAC3441643857628AE325418C299A9095AC3441F66C29B573E32541983F02087EAC34413080AA9F93E325419BAC3C0A86AC3441AC1F3B3DABE32541803A40B974AC3441E8CF3DB9C2E325413E3758C186AC3441D0AAB0E7F7E325410AAAA5429BAC3441BA971217DCE325418C721D6C98AC34415E2C5BB9D3E32541'); CREATE INDEX "bldg117862_geom_gist" ON "bldg117862" using gist ("geom" gist_geometry_ops); END;

0

Assuming that the orientation of the line segments is constant in a polygon, you could calculate the bearing (heading) of a vector perpendicular to each line segment. I don't have time right now to bash out to code but if you need the math, it can easily be supplied :-)

WolfOdrade
  • 2,748
  • 21
  • 24