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.
#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)
_______________________________________________
Proj mailing list
[hidden email]
http://lists.maptools.org/mailman/listinfo/proj
