Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

Michael Speth and i are in the process of fixing the PROJ.4 implementation of the HEALPix map projection (see http://trac.osgeo.org/proj/browser/trunk/proj/src/PJ_healpix.c for the current buggy implementation) and are running into some difficulties in the ellipsoidal case.

Here's the issue. We've implemented the spherical version of the HEALPix forward projection via

FORWARD(s_healpix_forward); /* sphere */
    return healpix_sphere(lp, P);
}

and that works fine. We now want to reuse healpix_sphere() in the ellipsoidal version of the HEALPix forward projection by simply replacing latitude lp.phi with authalic latitude and replacing the ellipsoid P with its authalic sphere. What's the best way to do this in PROJ.4?

We tried

FORWARD(e_healpix_forward); /* ellipsoid */
    lp.phi = auth_lat(lp.phi, P->e, 0);
    P->a = auth_radius(P->a, P->e);
    P->e = 0;
    return healpix_sphere(lp, P);
}

where auth_lat() and auth_radius() are our homemade authalic latitude and radius functions. But the numbers coming out aren't correct. We think that our attempt at changing P into its authalic sphere via the P->a and P->e assignments is failing. Any tips?

Michael also asked this question on the PROJ.4 forum (http://lists.maptools.org/pipermail/proj/2012-November/006465.html) last week but hasn't received a response.

Thanks for your attention.

share|improve this question

1 Answer 1

up vote 1 down vote accepted

Aha! I fixed the bug by setting P->a to authalic radius and P->ra to 1/P->a in the code below.

FORWARD(e_healpix_forward); // ellipsoid 
    lp.phi = auth_lat(P, lp.phi, 0);
    return healpix_sphere(lp);
}

...

ENTRY1(healpix, apa)
if (P->es) {
    P->apa = pj_authset(P->es); // For auth_lat().
    P->qp = pj_qsfn(1.0, P->e, P->one_es); // For auth_lat().
    P->a = P->a*sqrt(0.5*P->qp); // Set P->a to authalic radius.
    P->ra = 1.0/P->a;
    P->fwd = e_healpix_forward;
    P->inv = e_healpix_inverse; 
} else {
    P->fwd = s_healpix_forward;
    P->inv = s_healpix_inverse; 
}
ENDENTRY(P)

The introduction of P->apa is an inessential optimization for the auth_lat function that i use and was inspired by the cylindrical equal area projection source code PJ_cea.c.

For the full fixed PJ_healpix.c source code, which will soon be submitted as a patch to PROJ.4 (along with supporting test outputs), see https://www.dropbox.com/s/ocgu80fvkknhdev/PJ_healpix.c.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.