/******************************************************************************* NAME GCTP VERSION PROGRAMMER DATE ------- ---------- ---- T. Mittan 2-26-93 Conversion from FORTRAN to C S. Nelson 12-14-93 Added assignments to inunit and outunit for State Plane purposes. c.1.0 S. Nelson 9-15-94 Added outdatum parameter call. c.1.1 S. Nelson 11-94 Modified code so that UTM can accept any spheroid code. Changed State Plane legislated distance units, for NAD83, to be consistant with FORTRAN version of GCTP. Unit codes are specified by state laws as of 2/1/92. c.1.5 S. Nelson 1-98 Changed the name datum to spheroid. In initialize test, before the init functions were called for State Plane every time. Now, initialization will only be done if zone, spheroid, or the parameter array changes which is consistant with the other projections. For the UTM inverse projections the 1st 2 elements of a temporary projection array were being assigned to a lat and long that lies in that zone area. This has been eliminated. UTM will be initialized the same as the other projections for inverse transformations. For forward transformations the temporary array still exists. ALGORITHM REFERENCES 1. 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. 2. Snyder, John P. and Voxland, Philip M., "An Album of Map Projections", U.S. Geological Survey Professional Paper 1453 , United State Government Printing Office, Washington D.C., 1989. *******************************************************************************/ #include "cproj.h" #define TRUE 1 #define FALSE 0 static long iter = 0; /* First time flag */ static long inpj[MAXPROJ + 1]; /* input projection array */ static long indat[MAXPROJ + 1]; /* input dataum array */ static long inzn[MAXPROJ + 1]; /* input zone array */ static double pdin[MAXPROJ + 1][COEFCT];/* input projection parm array */ static long outpj[MAXPROJ + 1]; /* output projection array */ static long outdat[MAXPROJ + 1]; /* output dataum array */ static long outzn[MAXPROJ + 1]; /* output zone array */ static double pdout[MAXPROJ+1][COEFCT]; /* output projection parm array */ static long (*for_trans[MAXPROJ + 1])();/* forward function pointer array*/ static long (*inv_trans[MAXPROJ + 1])();/* inverse function pointer array*/ /* Table of unit codes as specified by state laws as of 2/1/92 for NAD 1983 State Plane projection, 1 = U.S. Survey Feet, 2 = Meters, 5 = International Feet */ static long NADUT[134] = {1, 5, 1, 1, 5, 1, 1, 1, 1, 2, 1, 1, 1, 2, 1, 1, 2, 2, 1, 1, 5, 2, 1, 2, 5, 1, 2, 2, 2, 1, 1, 1, 5, 2, 1, 5, 2, 2, 5, 2, 1, 1, 5, 2, 2, 1, 2, 1, 2, 2, 1, 2, 2, 2}; static long NAD83[134] = {101,102,5010,5300,201,202,203,301,302,401,402,403, 404,405,406,0000,501,502,503,600,700,901,902,903,1001,1002, 5101,5102,5103,5104,5105,1101,1102,1103,1201,1202,1301,1302, 1401,1402,1501,1502,1601,1602,1701,1702,1703,1801,1802,1900, 2001,2002,2101,2102,2103,2111,2112,2113,2201,2202,2203,2301, 2302,2401,2402,2403,2500,0000,0000,2600,0000,2701,2702,2703, 2800,2900,3001,3002,3003,3101,3102,3103,3104,3200,3301,3302, 3401,3402,3501,3502,3601,3602,3701,3702,3800,3900,0000,4001, 4002,4100,4201,4202,4203,4204,4205,4301,4302,4303,4400,4501, 4502,4601,4602,4701,4702,4801,4802,4803,4901,4902,4903,4904, 5001,5002,5003,5004,5005,5006,5007,5008,5009,5200,0000,5400}; void gctp(incoor,insys,inzone,inparm,inunit,inspheroid,ipr,efile,jpr,pfile, outcoor, outsys,outzone,outparm,outunit,outspheroid,fn27,fn83,iflg) double *incoor; /* input coordinates */ long *insys; /* input projection code */ long *inzone; /* input zone number */ double *inparm; /* input projection parameter array */ long *inunit; /* input units */ long *inspheroid; /* input spheroid */ long *ipr; /* printout flag for error messages. 0=screen, 1=file, 2=both*/ char *efile; /* error file name */ long *jpr; /* printout flag for projection parameters 0=screen, 1=file, 2 = both*/ char *pfile; /* error file name */ double *outcoor; /* output coordinates */ long *outsys; /* output projection code */ long *outzone; /* output zone */ double *outparm; /* output projection array */ long *outunit; /* output units */ long *outspheroid; /* output spheroid */ char fn27[]; /* file name of NAD 1927 parameter file */ char fn83[]; /* file name of NAD 1983 parameter file */ long *iflg; /* error flag */ { double x; /* x coordinate */ double y; /* y coordinate */ double factor; /* conversion factor */ double lon; /* longitude */ double lat; /* latitude */ double temp; /* temporary variable */ long i,j; /* loop counters */ long ininit_flag; /* input initilization flag */ long outinit_flag; /* output initilization flag */ long ind; /* temporary var used to find state plane zone */ long unit; /* temporary unit variable */ double temparr[COEFCT]; /* temporary projection array */ /* setup initilization flags and output message flags ---------------------------------------------------*/ ininit_flag = FALSE; outinit_flag = FALSE; *iflg = 0; *iflg = init(*ipr,*jpr,efile,pfile); if (*iflg != 0) return; /* check to see if initilization is required only the first 13 projection parameters are currently used. If more are added the loop should be increased. ---------------------------------------------------------*/ if (iter == 0) { for (i = 0; i < MAXPROJ + 1; i++) { inpj[i] = 0; indat[i] = 0; inzn[i] = 0; outpj[i] = 0; outdat[i] = 0; outzn[i] = 0; for (j = 0; j < COEFCT; j++) { pdin[i][j] = 0.0; pdout[i][j] = 0.0; } } ininit_flag = TRUE; outinit_flag = TRUE; iter = 1; } else { if (*insys != GEO) { if ((inzn[*insys] != *inzone) || (indat[*insys] != *inspheroid) || (inpj[*insys] != *insys)) { ininit_flag = TRUE; } else for (i = 0; i < 13; i++) if (pdin[*insys][i] != inparm[i]) { ininit_flag = TRUE; break; } } if (*outsys != GEO) { if ((outzn[*outsys] != *outzone) || (outdat[*outsys] != *outspheroid) || (outpj[*outsys] != *outsys)) { outinit_flag = TRUE; } else for (i = 0; i < 13; i++) if (pdout[*outsys][i] != outparm[i]) { outinit_flag = TRUE; break; } } } /* Check input and output projection numbers ------------------------------------------*/ if ((*insys < GEO) || (*insys > MAXPROJ)) { p_error("Insys is illegal","GCTP-INPUT"); *iflg = 1; return; } if ((*outsys < GEO) || (*outsys > MAXPROJ)) { p_error("Outsys is illegal","GCTP-OUTPUT"); *iflg = 2; return; } /* find the correct conversion factor for units ---------------------------------------------*/ unit = *inunit; /* use legislated unit table for State Plane -------------------------------------------*/ if ((*inspheroid == 0) && (*insys == SPCS) && (*inunit == STPLN_TABLE)) unit = FEET; if ((*inspheroid == 8) && (*insys == SPCS) && (*inunit == STPLN_TABLE)) unit = NADUT[(*inzone)/100]; /* find the factor unit conversions--all transformations are in radians or meters --------------------------------*/ if (*insys == GEO) *iflg = untfz(unit,RADIAN,&factor); else *iflg = untfz(unit,METER,&factor); if (*insys == SPCS) *inunit = unit; if (*iflg != 0) { close_file(); return; } x = incoor[0] * factor; y = incoor[1] * factor; /* Initialize inverse transformation ----------------------------------*/ if (ininit_flag) { inpj[*insys] = *insys; indat[*insys] = *inspheroid; inzn[*insys] = *inzone; for (i = 0;i < COEFCT; i++) pdin[*insys][i] = inparm[i]; /* Call the initialization function ----------------------------------*/ inv_init(*insys,*inzone,inparm,*inspheroid,fn27,fn83,iflg,inv_trans); if (*iflg != 0) { close_file(); return; } } /* Do actual transformations --------------------------*/ /* Inverse transformations ------------------------*/ if (*insys == GEO) { lon = x; lat = y; } else if ((*iflg = inv_trans[*insys](x, y, &lon, &lat)) != 0) { close_file(); return; } /* DATUM conversion should go here --------------------------------*/ /* The datum conversion facilities should go here */ /* Initialize forward transformation ----------------------------------*/ if (outinit_flag) { outpj[*outsys] = *outsys; outdat[*outsys] = *outspheroid; outzn[*outsys] = *outzone; for (i = 0;i < COEFCT; i++) pdout[*outsys][i] = outparm[i]; /* If the projection is UTM, copy to a temporary array. This way, the user does not have to enter the zone nor a lat/long within the zone. The program will calculate the zone from the input lat/long. ---------------------------------------------------------------------*/ if (*outsys == UTM) { for (i = 2; i < COEFCT; i++) temparr[i] = outparm[i]; temparr[0] = 0; temparr[1] = 0; if (outparm[0] == 0.0) { temparr[0] = pakr2dm(lon); temparr[1] = pakr2dm(lat); } else { temparr[0] = outparm[0]; temparr[1] = outparm[1]; } for_init(*outsys,*outzone,temparr,*outspheroid,fn27,fn83,iflg,for_trans); } else for_init(*outsys,*outzone,outparm,*outspheroid,fn27,fn83,iflg,for_trans); if (*iflg != 0) { close_file(); return; } } /* Forward transformations ------------------------*/ if (*outsys == GEO) { outcoor[0] = lon; outcoor[1] = lat; } else if ((*iflg = for_trans[*outsys](lon, lat, &outcoor[0], &outcoor[1])) != 0) { close_file(); return; } /* find the correct conversion factor for units ---------------------------------------------*/ unit = *outunit; /* use legislated unit table ----------------------------*/ if ((*outspheroid == 0) && (*outsys == SPCS) && (*outunit == STPLN_TABLE)) unit = 1; if ((*outspheroid == 8) && (*outsys == SPCS) && (*outunit == STPLN_TABLE)) unit = NADUT[(*outzone)/100]; if (*outsys == GEO) *iflg = untfz(RADIAN,unit,&factor); else *iflg = untfz(METER,unit,&factor); if (*outsys == SPCS) *outunit = unit; outcoor[0] *= factor; outcoor[1] *= factor; close_file(); return; }