¿Cómo arreglar PJ_healpix.c en PROJ.4?

Question

Michael Speth y yo estamos arreglando la implementación de PROJ.4 de la proyección cartográfica HEALPix (ver http://trac.osgeo.org/proj/browser/trunk/proj/src/PJ_healpix.c para la implementación actual con errores) y están encontrando algunas dificultades en el caso elipsoidal.

Esta es la cuestión. Hemos implementado la versión esférica de la proyección frontal de HEALPix mediante

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

y eso funciona bien. Ahora queremos reutilizar healpix_sphere() en la versión elipsoidal de la proyección frontal de HEALPix simplemente sustituyendo latitud lp.phi con la latitud auténtica y sustituyendo el elipsoide P con su esfera auténtica. ¿Cuál es la mejor manera de hacer esto en PROJ.4?

Intentamos

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);
}

donde auth_lat() y auth_radius() son nuestras funciones caseras de latitud y radio auténticos. Pero los números que salen no son correctos. Creemos que nuestro intento de cambiar P a su esfera auténtica a través de las asignaciones P->a y P->e está fallando. ¿Algún consejo?

Michael también hizo esta pregunta en el foro de PROJ.4 (http://lists.maptools.org/pipermail/proj/2012-November/006465.html) la semana pasada, pero no ha recibido respuesta.

Gracias por su atención.

Answer 1

Aha! He arreglado el error poniendo P->a como radio auténtico y P->ra como 1/P->a en el código de abajo.

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)

La introducción de P->apa es una optimización no esencial para la función auth_lat que utilizo y está inspirada en el código fuente de la proyección cilíndrica de área igual PJ_cea.c.

Para ver el código fuente completo de PJ_healpix.c, que pronto se enviará como parche a PROJ.4 (junto con los resultados de las pruebas), consulte https://www.dropbox.com/s/ocgu80fvkknhdev/PJ_healpix.c .