/* File: convert_albers.c Simple program to convert an x/y location, in meters, to a lat/lon location for NLCD 1992/2001 data that give the corner-point coordinates in x/y. The transformation assumes the x/y location is with respect to an Albers Equal Area projection using the NAD83 datum. When run, the program takes, as command-line arguments, the x and y locations; the corresponding latitude and longitude are printed to standard output. References: Snyder, John P., "Map Projections -- A Working Manual", USGS Professional Paper 1395, United State Government Printing Office, Washington D.C., 1987. Michael G. Duda, NCAR/MMM */ #include #include #include #define RE_NAD83 6378137.0 #define E_NAD83 0.0818187034 /* Eccentricity */ #define D2R 0.017453293 /* PI/180 */ int main(int argc, char ** argv) { double x, y, lat, lon; double truelat1, truelat2, stand_lon, lat1, lon1; double rho, rho0, theta, m1, m2, sinphi, C, nc, q1, q2, q0, q, beta, a, b, c; /* True latitudes and standard longitude used for NLCD */ truelat1 = 29.5; truelat2 = 45.5; stand_lon = -96.0; /* Lat and lon of origin */ lat1 = 23.0; lon1 = -96.0; if (argc != 3) { fprintf(stderr,"Usage: convert_albers x y\n\n"); return 1; } x = atof(argv[1]); y = atof(argv[2]); m1 = cos(truelat1*D2R)/sqrt(1.0-pow((E_NAD83*sin(truelat1*D2R)),2.0)); m2 = cos(truelat2*D2R)/sqrt(1.0-pow((E_NAD83*sin(truelat2*D2R)),2.0)); sinphi = sin(truelat1*D2R); q1 = (1.0-pow(E_NAD83,2.0))*((sinphi/(1.0-pow((E_NAD83*sinphi),2.0))) - 1.0/(2.0*E_NAD83) * log((1.0-E_NAD83*sinphi)/(1.0+E_NAD83*sinphi))); sinphi = sin(truelat2*D2R); q2 = (1.0-pow(E_NAD83,2.0))*((sinphi/(1.0-pow((E_NAD83*sinphi),2.0))) - 1.0/(2.0*E_NAD83) * log((1.0-E_NAD83*sinphi)/(1.0+E_NAD83*sinphi))); nc = (pow(m1,2.0) - pow(m2,2.0)) / (q2 - q1); C = pow(m1,2.0) + nc*q1; sinphi = sin(lat1*D2R); q0 = (1.0-pow(E_NAD83,2.0))*((sinphi/(1.0-pow((E_NAD83*sinphi),2.0))) - 1.0/(2.0*E_NAD83) * log((1.0-E_NAD83*sinphi)/(1.0+E_NAD83*sinphi))); rho0 = RE_NAD83 * sqrt(C - nc * q0) / nc; rho = sqrt(pow(x,2.0) + pow((rho0 - y),2.0)); q = (C - pow((rho*nc/RE_NAD83),2.0)) / nc; beta = asin(q/(1.0 - log((1.0-E_NAD83)/(1.0+E_NAD83))*(1.0-pow(E_NAD83,2.0))/(2.0*E_NAD83))); a = 1./3.*pow(E_NAD83,2.) + 31./180.*pow(E_NAD83,4.) + 517./5040.*pow(E_NAD83,6.); b = 23./360.*pow(E_NAD83,4.) + 251./3780.*pow(E_NAD83,6.); c = 761./45360.*pow(E_NAD83,6.); theta = atan2(x, (rho0 - y)); lat = (beta + a*sin(2.0*beta) + b*sin(4.0*beta) + c*sin(6.0*beta))/D2R; lon = stand_lon + (theta/D2R)/nc; printf("known_lat = %f\n",lat); printf("known_lon = %f\n",lon); return 0; }