[GRASS-dev] /lib/gis/area_poly1.c: algorithm poorly described, new comments for consideration please

[Richard Roger]

Dear Dev list,
 
I hope this is the correct place to post a request for an update to two comments in /lib/gis/area_poly1.c, and replace it with something more informative.
 
My remarks are based on reading /lib/gis/area_poly1.c from     grass-7.0.svn_src_snapshot_2013_07_13.
 
This file has routines for computing areas of polygons on the ellipsoid.  The algorithm is hardly described, and it is not trivial to work out what is going on.  My comment addresses that deficiency.
 
Regards,
Richard Roger
 
 
1.     The comment before the routine G_begin_ellipsoid_polygon_area reads:
 
 * \brief Begin area calculations.
 *
 * This initializes the polygon area calculations for the
 * ellipsoid with semi-major axis <i>a</i> (in meters) and ellipsoid
 * eccentricity squared <i>e2</i>.
 *
 * \param a semi-major axis
 * \param e2 ellipsoid eccentricity
 
For the comment to be self-consistent, the last line should be:
 
 * \param e2 ellipsoid eccentricity squared
 
2.    The comment before the routine G_ellipsoid_polygon_area describes briefly what it
computes:
 
 * <b>Note:</b> This routine assumes grid lines on the connecting the 
 * vertices (as opposed to geodesics).
 
I find this comment less than useful.  The method this routine implements for computing area is
not one that I have been able to find explicitly in the literature, i.e., I did not recognize the Qbar
function nor the QBarA, QbarB, etc., constants.  As such, it is not at all obvious just what area this
routine computes.  With some work, I have been able to figure out what is going on.  I suggest
the Note be replaced with, e.g.,
 
 * <b>Note:</b> This routine computes the area of a polygon on the
 * ellipsoid.  The sides of the polygon are, in general, not geodesics.
 * Each side is actually defined by a linear relationship between
 * latitude and longitude, i.e., on a rectangular/equidistant
 * cylindrical/Plate Carr{'e}e grid, the side would appear as a
 * straight line.  For two consecutive vertices of the polygon, 
 * (lat_1, long1) and (lat_2,long_2), the line joining them (i.e., the
 * polygon's side) is defined by:
 *                                     lat_2  -  lat_1 
 *     lat = lat_1 + (long - long_1) * ---------------
 *                                     long_2 - long_1
 * where long_1 < long < long_2.  This is not a straight line on 
 * the ellipsoid.
 *   The values of QbarA, etc., are determined by the integration of
 * the Q function.  Into www.integral-calculator.com, paste this 
 * expression :
 *
 * sin(x)+ (2/3)e^2(sin(x))^3 + (3/5)e^4(sin(x))^5 + (4/7)e^6(sin(x))^7
 *
 * and you'll get their values.  (Last checked 30 Oct 2013). 
 *
 * This function correctly computes (within the limits of the series
 * approximation) the area of a quadrilateral on the ellipsoid when
 * two of its sides run along meridians and the other two sides run
 * along parallels of latitude.
 
 
 
 


Dr. R. E. Roger| Senior Spatial Analyst
NSW Department of Primary Industries |NSW Office of Water
181 - 183 Anson Street| Orange NSW 2800 | PO Box 53 | Orange NSW 2800
T: 02 6363 8729 |F: 02 6361 3839 |E: Richard.Roger at water.nsw.gov.au
W: www.water.nsw.gov.au ( http://www.trade.nsw.gov.au/ )


-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/grass-dev/attachments/20131031/a99567e3/attachment.html>