home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 3 / RISC_DISC_3.iso / resources / etexts / gems / gemsi / albers.c < prev    next >
C/C++ Source or Header  |  1992-09-15  |  4KB  |  161 lines
  1. /*
  2. Albers Equal-Area Conic Map Projection
  3. by Paul Bame
  4. from "Graphics Gems", Academic Press, 1990
  5. */
  6.  
  7. /*
  8.  * Albers Conic Equal-Area Projection
  9.  * Formulae taken from "Map Projections Used by the U.S. 
  10.  * Geological Survey" Bulletin #1532
  11.  *
  12.  * Equation reference numbers and some variable names taken
  13.  * from the reference.
  14.  * To use, call albers setup() once and then albers_project()
  15.  * for each coordinate pair of interest.
  16. */
  17.  
  18. #include "GraphicsGems.h"        
  19. #include <stdio.h>
  20. #include <math.h>
  21.  
  22. /*
  23.  * This is the Clarke 1866 Earth spheroid data which is often
  24.  * used by the USGS - there are other spheroids however - see the
  25.  * book.
  26.  */
  27.  
  28. /*
  29.  * Earth radii in different units */
  30. #define CONST_EradiusKM (6378.2064)     /* Kilometers */
  31. #define CONST_EradiusMI (CONST_EradiusKM/1.609)    /* Miles */
  32. #define CONST_Ec        (0.082271854)    /* Eccentricity */
  33. #define CONST_Ecsq    (0.006768658)    /* Eccentricity squared */
  34.  
  35.  
  36.  
  37. /*
  38.  * To keep things simple, assume Earth radius is 1.0. Projected
  39.  * coordinates (X and Y obtained from albers project ()) are
  40.  *  dimensionless and may be multiplied by any desired radius
  41.  *  to convert to desired units (usually Kilometers or Miles).
  42. */
  43. #define CONST_Eradius 1.0
  44.  
  45. /* Pre-computed variables */
  46. static double middlelon;        /* longitude at center of map */
  47. static double bigC, cone_const, r0;    /* See the reference */
  48.  
  49. static
  50. calc_q_msq(lat, qp, msqp)
  51. double lat;
  52. double *qp;
  53. double *msqp;
  54. /*
  55.  * Given latitude, calculate 'q' [eq 3-12] 
  56.  * if msqp is != NULL, m^2  [eq. 12-15].
  57. */
  58. {
  59.     double s, c, es;
  60.  
  61.     s = sin(lat);
  62.     es = s * CONST_Ec;
  63.  
  64.     *qp = (1.0 - CONST_Ecsq) * ((s / (1 - es * es))-
  65.         (1 / (2 * CONST_Ec)) * log((1 - es) / (1 + es)));
  66.  
  67.     if (msqp != NULL)
  68.     {
  69.         c = cos(lat);
  70.         *msqp = c * c/ (1 - es * es);
  71.     }
  72. }
  73.  
  74.  
  75.  
  76.  
  77. albers_setup(southlat, northlat, originlon, originlat)
  78. double southlat, northlat;
  79. double originlon;
  80. double originlat;
  81. /*
  82.  * Pre-compute a bunch of variables which are used by the 
  83.  * albers_project()
  84.  *
  85.  * southlat     Southern latitude for Albers projection
  86.  * northlat    Northern latitude for Albers projection
  87.  * originlon    Longitude for origin of projected map
  88.  * originlat    Latitude for origin of projected map -
  89.  *                often (northlat + southlat) / 2
  90. */
  91. {
  92.     double q1, q2, q0;
  93.     double m1sq, m2sq;
  94.  
  95.     middlelon = originlon;
  96.     
  97.     calc_q_msq(southlat, &q1, &m1sq);
  98.     calc_q_msq(northlat, &q2, &m2sq);
  99.     calc_q_msq(originlat, &q0, (double *)NULL);
  100.  
  101.     cone_const = (m1sq - m2sq) / (q2 - q1);
  102.     bigC = m1sq + cone_const * q1;
  103.     r0 = CONST_Eradius * sqrt(bigC - cone_const *q0) / cone_const;
  104. }
  105.  
  106. /***************************************************************/
  107.  
  108. albers_project(lon, lat, xp, yp)
  109. double lon, lat;
  110. double *xp, *yp;
  111. /*
  112.  * Project lon/lat (in radians) according to albers_setup and 
  113.  * return the results via xp, yp. Units of *xp and *yp are same
  114.  * as the units used by CONST_Eradius.
  115. */
  116. {
  117.     double q, r, theta;
  118.     calc_q_msq(lat, &q, (double *)NULL);
  119.     theta = cone_const * (lon -middlelon);
  120.     r = CONST_Eradius * sqrt(bigC - cone_const * q) / cone_const;
  121.     *xp = r * sin(theta);
  122.     *yp = r0 - r * cos(theta);
  123. }
  124.  
  125. #ifdef TESTPROGRAM
  126.  
  127. /*
  128.  * Test value from the USGS book. Because of roundoff, the 
  129.  * actual values are printed for visual inspection rather 
  130.  * than guessing what constitutes "close enough".
  131. */
  132. /* Convert a degrees, minutes pair to radians */
  133. #define DEG_MIN2RAD(degrees, minutes) \
  134.     ((double) ((degrees + minutes / 60.0) * M_PI / 180.0))
  135.  
  136. #define Nlat DEG_MIN2RAD(29, 30)    /* 29 degrees, 30' North Lat */
  137. #define Slat DEG_MIN2RAD(45, 30)    
  138. #define Originlat DEG_MIN2RAD(23, 0)    
  139. #define Originlon DEG_MIN2RAD(-96, 0) /* West longitude is negative */
  140.  
  141. #define Testlat DEG_MIN2RAD(35, 0)
  142. #define Testlon DEG_MIN2RAD(-75, 0)
  143.  
  144. #define TestX 1885.4727
  145. #define TestY 1535.9250
  146.  
  147. main()
  148. {
  149.     int i;
  150.     double x, y;
  151.  
  152. /* Setup is also from USGS book test set */
  153.     albers_setup(Slat, Nlat, Originlon, Originlat);
  154.  
  155.     albers_project(Testlon, Testlat, &x, &y);
  156.     printf("%.41f, %.41f =?= %.41f, %.41f/n",
  157.         x * CONST_EradiusKM, y * CONST_EradiusKM,
  158.         TestX, TestY);
  159. }
  160. #endif        /* TESTPROGRAM */
  161.