/* This routine maps Latitude/Longitude to Line sample coordinates in the Conterminous U.S. AVHRR data set on CD-ROM. This data set is in the Lambert Azimuthal Equal Area projection. This function 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 Latitude,Longitude to Line,Sample\n\n"); /* Initialize the Lambert Azimuthal Equal Area projection ------------------------------------------------------*/ lamaz_init(radius, center_lon, center_lat); /* Process each point ------------------*/ printf("\nAt the > prompt, enter Latitude Longitude in decimal degrees\n"); printf("Separate numbers with blanks. Terminate with -999 -999\n"); printf("Results returned are Line,Sample coordinates.\n\n"); for (;;) { printf("> "); scanf("%lf %lf", &in1, &in2); if (in1 == -999.0) break; lat = in1 * D2R; lon = in2 * D2R; lamaz_forward(lon, lat, &x, &y); /* Here, x and y are projection coordinates in the Lambert Azimuthal Equal Area projection. They are then converted to line, sample coordinates. For this dataset, the center of the pixel (1,1) is (-2050000,75200) and the pixel size is 1000.0 meters. --------------------------------*/ line = (752000.0 - y) / 1000.0 + 1.0; samp = (x + 2050000.0) / 1000.0 + 1.0; printf("%15.3f %15.3f\n",line,samp); } } /* 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 forward equations--mapping lat,long to x,y -----------------------------------------------------------------------*/ lamaz_forward(lon, lat, x, y) double lon; /* (I) Longitude */ double lat; /* (I) Latitude */ double *x; /* (O) X projection coordinate */ double *y; /* (O) Y projection coordinate */ { double adjust_lon(); /* Function to adjust longitude to -180 - 180 */ double delta_lon; /* Delta longitude (Given longitude - center */ double sin_delta_lon; /* Sine of the delta longitude */ double cos_delta_lon; /* Cosine of the delta longitude */ double sin_lat; /* Sine of the given latitude */ double cos_lat; /* Cosine of the given latitude */ double g; double ksp; /* Forward equations -----------------*/ delta_lon = adjust_lon(lon - lon_center); sincos(lat, &sin_lat, &cos_lat); sincos(delta_lon, &sin_delta_lon, &cos_delta_lon); g = sin_lat_o * sin_lat + cos_lat_o * cos_lat * cos_delta_lon; if (g == -1.0) { printf("Point projects to a circle of radius = %lf\n", 2.0 * R); return; } ksp = R * sqrt(2.0 / (1.0 + g)); *x = ksp * cos_lat * sin_delta_lon; *y = ksp * (cos_lat_o * sin_lat - sin_lat_o * cos_lat * cos_delta_lon); 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)