[Proj] Need help converting Goode to latlon

[David Turover]

On Thu, Jul 23, 2009 at 07:37:55PM -0400, Gerald I. Evenden wrote:
> If the interupted version of Goode is used then *NO* standard version of 
> [l]proj will handle it without the user segmenting the data into the 
> subregions of the interuptions.

Thanks. If attachments work on this list, this post should carry a
rough modified version of PJ_goode which seems to work for at least
the three regions that I'm working in. I have not tested it more fully. 

 - David T. 
-------------- next part --------------
#ifndef lint
static const char SCCSID[]="@(#)PJ_goode.c	4.1 94/02/15     GIE     REL";
#endif
#define PROJ_PARMS__ \
	struct PJconsts	*sinu; \
	struct PJconsts	*moll;
#define PJ_LIB__
#include	<projects.h>
PROJ_HEAD(goode_interrupted, "Goode Homolosine Interrupted") "\n\tPCyl, Sph.";
	C_NAMESPACE PJ
*pj_sinu(PJ *), *pj_moll(PJ *);
#define Y_COR		0.05280
#define PHI_LIM	.71093078197902358062

/* Goode interruption parameters from:
 * http://edc2.usgs.gov/1KM/goodesarticle.php
 * A goode map for visual reference is Figure 5 from here: 
 * http://proceedings.esri.com/library/userconf/proc98/proceed/to850/pap844/p844.htm
 */
static const double goode_center_x_for_region[] = {
	 0, /* Default. */
	/* Use radians. */
	PI * -100.0 / 180.0,  /* Region 1 */
	PI *   30.0 / 180.0,  /* Region 2 */
	PI * -100.0 / 180.0,  /* Region 3 */
	PI *   30.0 / 180.0,  /* Region 4 */
	PI * -160.0 / 180.0,  /* Region 5 */
	PI * - 60.0 / 180.0,  /* Region 6 */
	PI *   20.0 / 180.0,  /* Region 7 */
	PI *  140.0 / 180.0,  /* Region 8 */
	PI * -160.0 / 180.0,  /* Region 9 */
	PI * - 60.0 / 180.0,  /* Region 10*/
	PI *   20.0 / 180.0,  /* Region 11*/
	PI *  140.0 / 180.0,  /* Region 12*/
};

int goode_get_region_from_lat_lon(double lat, double lon){
	if(lon > 0 ){ /* North */
		if(lon > 0.45259259259 * PI){ /* Norther */
			if(lat < -.2222222222 * PI){ 
				return 1; 
			} else {
				return 2;
			}
		} else {
			/* Note: Duped code */
			if(lat < -.2222222222 * PI){ 
				return 3; 
			} else {
				return 4;
			}
		}

	} else { /* South */
		if(lon > -0.45259259259 * PI){ /* First section south of equator */ 
			if(lat < .1111111111 * PI){ 
				if(lat < -.8888888888 * PI){ 
					return 5;
				} else {
					return 6;
				}
			} else {
				if(lat < .7777777777 * PI){ 
					return 7;
				} else {
					return 8;
				}
			}
		} else {
			/* Note: Duped code */
			if(lat < .1111111111 * PI){ 
				if(lat < -.8888888888 * PI){ 
					return 9;
				} else {
					return 10;
				}
			} else {
				if(lat < .7777777777 * PI){ 
					return 11;
				} else {
					return 12;
				}
			}
		}
	}
	return 0;
}

/* Macro alternative to the function. */
/* Note: Using a macro does not appear to save much time. */
#define GOODE_REGION_FROM_LATLON(lat, lon) \
( \
	(lon > 0.0) \
	?( \
		1 + \
		( 0x01 & (lat > (- PI * 40.0 / 180.0)) \
		| 0x02 & (lon < (PI * (40.0 + 44.0/60.0) / 90.0 )) \
		) \
	):( \
		5 +  \
		( 0x01 & (lat < (PI * 140.0 / 180.0)) \
		| 0x02 & (lat > (PI * 20.0 / 180.0)) \
		| 0x04 & (lon < (- PI * (40.0 + 44.0/60.0) / 90.0 )) \
		)\
	)\
)
		

FORWARD(s_forward); /* spheroid */
	/* int goode_current_region = goode_get_region_from_lat_lon(lp.lam, lp.phi); */
	int goode_current_region = GOODE_REGION_FROM_LATLON(lp.lam, lp.phi);
	
	fprintf(stderr, "In forward, sinu %x moll %x reg %d, lam0 %.4f, l,p are: %.4f, %.4f\n",
	P->sinu, P->moll, 
	goode_current_region, P->lam0, lp.lam, lp.phi);

	lp.lam -= goode_center_x_for_region[goode_current_region];

	if (fabs(lp.phi) <= PHI_LIM)
		xy = P->sinu->fwd(lp, P->sinu);
	else {
		xy = P->moll->fwd(lp, P->moll);
		xy.y -= lp.phi >= 0.0 ? Y_COR : -Y_COR;
	}

	fprintf(stderr, "In forward after mods lam0 %.4f, x,y are: %.4f, %.4f lp is %.4f, %.4f\n", 
	P->lam0, xy.x, xy.y, lp.lam, lp.phi);

	xy.x += goode_center_x_for_region[goode_current_region];
	return (xy);
}
INVERSE(s_inverse); /* spheroid */
	/* int goode_current_region = goode_get_region_from_lat_lon(xy.x, xy.y);*/
	int goode_current_region = GOODE_REGION_FROM_LATLON(xy.x, xy.y);

	fprintf(stderr, "In inverse, reg %d, x,y are: %.4f, %.4f\n",
	goode_current_region, xy.x, xy.y);

	xy.x -= goode_center_x_for_region[goode_current_region];

	if (fabs(xy.y) <= PHI_LIM)
		lp = P->sinu->inv(xy, P->sinu);
	else {
		xy.y += xy.y >= 0.0 ? Y_COR : -Y_COR;
		lp = P->moll->inv(xy, P->moll);
	}

	fprintf(stderr, "In inverse after mods lam0 %.4f, x,y are: %.4f, %.4f lp is %.4f, %.4f\n", 
	P->lam0, xy.x, xy.y, lp.lam, lp.phi);

	lp.lam += goode_center_x_for_region[goode_current_region];
	return (lp);
}
FREEUP;
	if (P) {
		if (P->sinu)
			(*(P->sinu->pfree))(P->sinu);
		if (P->moll)
			(*(P->moll->pfree))(P->moll);
		pj_dalloc(P);
	}
}
ENTRY2(goode_interrupted, sinu, moll)
	P->es = 0.;
	if (!(P->sinu = pj_sinu(0)) || !(P->moll = pj_moll(0)))
		E_ERROR_0;
	if (!(P->sinu = pj_sinu(P->sinu)) || !(P->moll = pj_moll(P->moll)))
		E_ERROR_0;
	P->fwd = s_forward;
	P->inv = s_inverse;
ENDENTRY(P)