/******************************************************************************* 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. 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 "gctpc/cproj.h" #include "gctpc/proj.h" #include "gctpc/report.h" #include "gctpc/untfz.h" #include "gctpc/gctp.h" #include "gctpc/for_init.h" #include "gctpc/inv_init.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][15]; /* 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][15]; /* output projection parm array */ static ForwardTransFunc for_trans[MAXPROJ + 1];/* forward function pointer array*/ static InverseTransFunc 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(double *incoor, long int *insys, long int *inzone, double *inparm, long int *inunit, long int *indatum, long int *ipr, char *efile, long int *jpr, char *pfile, double *outcoor, long int *outsys, long int *outzone, double *outparm, long int *outunit, long int *outdatum, char *fn27, char *fn83, long int *iflg) /* input coordinates */ /* input projection code */ /* input zone number */ /* input projection parameter array */ /* input units */ /* input datum */ /* printout flag for error messages. 0=screen, 1=file, 2=both*/ /* error file name */ /* printout flag for projection parameters 0=screen, 1=file, 2 = both*/ /* error file name */ /* output coordinates */ /* output datum */ /* output projection code */ /* output zone */ /* output projection array */ /* output units */ /* file name of NAD 1927 parameter file */ /* file name of NAD 1983 parameter file */ /* error flag */ { double x; /* x coordinate */ double y; /* y coordinate */ double factor; /* conversion factor */ double lon; /* longitude */ double lat; /* latitude */ double temp; /* dummy 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 dummy[15]; /* 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 < 15; 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] != *indatum) || (inpj[*insys] != *insys) || (*insys == 2)) { 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] != *outdatum) || (outpj[*outsys] != *outsys) || (*outsys == 2)) { 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 < 0) || (*insys > MAXPROJ)) { p_error("Insys is illegal","GCTP-INPUT"); *iflg = 1; return; } if ((*outsys < 0) || (*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 ((*indatum == 0) && (*insys == 2) && (*inunit == 6)) unit = 1; if ((*indatum == 8) && (*insys == 2) && (*inunit == 6)) unit = NADUT[(*inzone)/100]; /* find the factor unit conversions --------------------------------*/ if (*insys == GEO) *iflg = untfz(unit,0,&factor); else *iflg = untfz(unit,2,&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] = *indatum; inzn[*insys] = *inzone; for (i = 0;i < 15; i++) pdin[*insys][i] = inparm[i]; if (*insys == 1) { for( i = 2; i < 15; i++) dummy[i] = inparm[i]; dummy[0] = 0; dummy[1] = 0; if ((*inzone != 0) || (inparm[0] == 0.0)) { dummy[0] = 1.0e6 * (double)(6 * *inzone -183); if ( *inzone >= 0) dummy[1] = 4.0e7; else dummy[1] = -4.0e7; } else { dummy[0] = inparm[0]; dummy[1] = inparm[1]; } inv_init(*insys,*inzone,dummy,*indatum,fn27,fn83,iflg,inv_trans); } else inv_init(*insys,*inzone,inparm,*indatum,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] = *outdatum; outzn[*outsys] = *outzone; for (i = 0;i < 15; i++) pdout[*outsys][i] = outparm[i]; if (*outsys == 1) { for (i = 2; i < 15; i++) dummy[i] = outparm[i]; dummy[0] = 0; dummy[1] = 0; if (outparm[0] == 0.0) { dummy[0] = pakr2dm(lon); dummy[1] = pakr2dm(lat); } else { dummy[0] = outparm[0]; dummy[1] = outparm[1]; } for_init(*outsys,*outzone,dummy,*outdatum,fn27,fn83,iflg,for_trans); } else for_init(*outsys,*outzone,outparm,*outdatum,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 ((*outdatum == 0) && (*outsys == 2) && (*outunit == 6)) unit = 1; if ((*outdatum == 8) && (*outsys == 2) && (*outunit == 6)) unit = NADUT[(*outzone)/100]; if (*outsys == GEO) *iflg = untfz(0,unit,&factor); else *iflg = untfz(2,unit,&factor); if (*outsys == SPCS) *outunit = unit; outcoor[0] *= factor; outcoor[1] *= factor; close_file(); return; }