/* This routine maps Line/Sample coordinates in the Conterminous U.S. AVHRR data set on CD-ROM 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: Line Sample Latitude Longitude --------------- ---------------------- 1 1 48.400507 -128.521181 1 4587 46.702829 -65.403539 2889 1 23.589215 -119.968362 2889 4587 22.485019 -75.420050 .5 .5 48.403055 -128.530059 .5 4587.5 46.704899 -65.394649 2889.5 .5 23.583758 -119.972290 2889.5 4587.5 22.479392 -75.416353 -999 -999 Please note that the second group of four test points returns the Latitude, Longitude values of the image minimum bounding rectangle. This routine bases it's calculations on pixel-centered coordinates, whereas the minimum bounding rectangle in the README file gives the coordinates of the pixel corners. 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; double line, samp; printf("Routine to map Line/Sample coordinates to Latitude/Longitude.\n\n"); /* Initialize the Lambert Azimuthal Equal Area projection ------------------------------------------------------*/ lamaz_init(radius, center_lon, center_lat); /* Process each point ------------------*/ printf("\nAt the > prompt, enter Line Sample coordinates.\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", &line, &samp); if (line == -999.0) break; /* Convert Line/Sample to projection space meters and convert to lat/long. The center of the pixel is assumed. -----------------------------------*/ in1 = -2050000.0 + ((samp - 1) * 1000.0); in2 = 752000.0 - ((line - 1) * 1000.0); 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)