/* This routine maps Lambert Azimuthal Equal Area projection coordinates to Latitude/Longitude. It was complied and run using the C compiler on SunOS 4.1 (UNIX) and on MS-DOS 3.x using TURBO C++. Results on these machines were accurate at the time of testing, but they are by no means guaranteed! Some test data: Latitude/Longitude Projection X,Y ------------------ -------------- 20.100000 -129.900000 -3127595.117 -2222461.848 20.100000 -119.900000 -2107546.588 -2514593.914 20.100000 -109.900000 -1056159.350 -2689483.336 20.100000 -99.900000 10693.712 -2747009.429 20.100000 -89.900000 1077390.960 -2687135.457 20.100000 -79.900000 2128308.259 -2509899.685 20.100000 -69.900000 3147564.983 -2215424.707 20.100000 -59.900000 4118744.741 -1803950.241 30.100000 -129.900000 -2830213.099 -1153081.051 30.100000 -119.900000 -1909925.977 -1430783.486 30.100000 -109.900000 -957941.342 -1597312.839 30.100000 -99.900000 9701.938 -1652133.370 30.100000 -89.900000 977187.412 -1595075.846 30.100000 -79.900000 1928696.405 -1426316.480 30.100000 -69.900000 2848183.241 -1146398.866 30.100000 -59.900000 3719136.845 -756304.406 40.100000 -129.900000 -2476574.541 -87188.349 40.100000 -119.900000 -1673888.753 -341510.104 40.100000 -109.900000 -840327.554 -494327.530 40.100000 -99.900000 8513.313 -544683.714 40.100000 -89.900000 857200.107 -492273.226 40.100000 -79.900000 1690297.453 -337414.058 40.100000 -69.900000 2492204.689 -81076.846 40.100000 -59.900000 3246990.758 275024.244 50.100000 -129.900000 -2070732.002 966649.168 50.100000 -119.900000 -1401942.702 744697.524 50.100000 -109.900000 -704507.500 611013.737 50.100000 -99.900000 7139.659 566911.169 50.100000 -89.900000 718643.394 612812.385 50.100000 -79.900000 1415647.521 748277.476 50.100000 -69.900000 2083715.752 971974.569 50.100000 -59.900000 2708260.488 1281590.114 -999.000000 -999.000000 -999.000 -999.000 by DRS/TGS/ISS/EROS Data Center, May 1991 Algorithm references: 1. "New Equal-Area Map Projections for Noncircular Regions", John P. Snyder, The American Cartographer, Vol 15, No. 4, October 1988, pp. 341-355. 2. Snyder, John P., "Map Projections--A Working Manual", U.S. Geological Survey Professional Paper 1395 (Supersedes USGS Bulletin 1532), United State Government Printing Office, Washington D.C., 1987. 3. "Software Documentation for GCTP General Cartographic Transformation Package", U.S. Geological Survey National Mapping Division, May 1982. -------------------------------------------------------------------------*/ #include #include #define PI 3.1415926536 #define HALF_PI 1.5707963268 #define TWO_PI 6.2831853072 #define EPSLN 1.0e-10 #define R2D 57.2957795131 #define D2R 0.01745329252 /* Variables common to all subroutines in this code file -----------------------------------------------------*/ static double lon_center; /* Center longitude (projection center) */ static double lat_center; /* Center latitude (projection center) */ static double R; /* Radius of the earth (sphere) */ static double sin_lat_o; /* Sine of the center latitude */ static double cos_lat_o; /* Cosine of the center latitude */ void main() { double lat,lon; double x,y; double radius = 6370997.0; /* Radius of sphere in meters */ double center_lon = -1.74532925199; /* Center long is -100.0 degrees */ double center_lat = 0.785398163; /* Center latitude is 45.0 degrees */ double in1, in2; printf("Routine to map Lambert Azimuthal EA projection X,Y to Lat/Long\n\n"); /* Initialize the Lambert Azimuthal Equal Area projection ------------------------------------------------------*/ lamaz_init(radius, center_lon, center_lat); /* Process each point ------------------*/ printf("\nAt the > prompt, enter Projection X Y in meters \n"); printf("Separate numbers with blanks. Terminate with -999 -999\n"); printf("Results are returned Latitude Longitude in decimal degrees\n\n"); for (;;) { printf("> "); scanf("%lf %lf", &in1, &in2); if (in1 == -999.0) break; lamaz_inverse(in1, in2, &lon, &lat); lat *= R2D; lon *= R2D; printf("%12.6f %12.6f\n",lat,lon); } } /* Initialize the Lambert Azimuthal Equal Area projection ------------------------------------------------------*/ lamaz_init(r, center_long, center_lat) double r; /* (I) Radius of the earth (sphere) */ double center_long; /* (I) Center longitude */ double center_lat; /* (I) Center latitude */ { /* Place parameters in static storage for common use -------------------------------------------------*/ R = r; lon_center = center_long; lat_center = center_lat; sincos(center_lat, &sin_lat_o, &cos_lat_o); /* Report parameters to the user -----------------------------*/ printf("Lambert Azimuthal Equal-Area Projection Parameters:\n"); printf("---------------------------------------------------\n"); printf(" Radius of Sphere: %lf meters\n", r); printf(" Longitude of Center: %lf degrees\n", center_long * R2D); printf(" Latitude of Center: %lf degrees\n", center_lat * R2D); return; } /* Lambert Azimuthal Equal Area inverse equations--mapping x,y to lat,long -----------------------------------------------------------------------*/ lamaz_inverse(x, y, lon, lat) double x; /* (I) X projection coordinate */ double y; /* (I) Y projection coordinate */ double *lon; /* (O) Longitude */ double *lat; /* (O) Latitude */ { double adjust_lon(); /* Function to adjust longitude to -180 - 180 */ double Rh; double z; /* Great circle dist from proj center to given point */ double sin_z; /* Sine of z */ double cos_z; /* Cosine of z */ double temp; /* Re-used temporary variable */ /* Inverse equations -----------------*/ Rh = sqrt(x * x + y * y); temp = Rh / (2.0 * R); if (temp > 1) { printf("Input data error\n"); return; } z = 2.0 * asin(temp); sincos(z, &sin_z, &cos_z); *lon = lon_center; if (fabs(Rh) > EPSLN) { *lat = asin(sin_lat_o * cos_z + cos_lat_o * sin_z * y / Rh); temp = fabs(lat_center) - HALF_PI; if (fabs(temp) > EPSLN) { temp = cos_z - sin_lat_o * sin(*lat); if(temp!=0.0)*lon=adjust_lon(lon_center+atan2(x*sin_z*cos_lat_o,temp*Rh)); } else if (lat_center < 0.0) *lon = adjust_lon(lon_center - atan2(-x, y)); else *lon = adjust_lon(lon_center + atan2(x, -y)); } else *lat = lat_center; return; } /* Function to calculate the sine and cosine in one call. Some computer systems have implemented this function, resulting in a faster implementation than calling each function separately. It is provided here for those computer systems which don`t implement this function ----------------------------------------------------*/ #ifndef sun sincos(val, sin_val, cos_val) double val; double *sin_val; double *cos_val; { *sin_val = sin(val); *cos_val = cos(val); return; } #endif /* Function to return the sign of an argument ------------------------------------------*/ sign(x) double x; { if (x < 0.0) return(-1); else return(1);} /* Function to adjust longitude to -180 to 180 -------------------------------------------*/ double adjust_lon(x) double x;{x=(fabs(x)