PJ_aea.cs

[Pat28]

I just noticed a multi threading issue in AlbersEqualArea. In both FORWARD and INVERSE the member P->rho is modified. Meaning that if you happened to use a multi-threaded app that happened to use the same instance, you get unpredictable results. My guess is that P->rho should be a local variable.

FORWARD(e_forward); /* ellipsoid & spheroid */
    if ((P->rho = P->c - (P->ellips ? P->n * pj_qsfn(sin(lp.phi),
        P->e, P->one_es) : P->n2 * sin(lp.phi))) < 0.) F_ERROR
    P->rho = P->dd * sqrt(P->rho);
    xy.x = P->rho * sin( lp.lam *= P->n );
    xy.y = P->rho0 - P->rho * cos(lp.lam);
    return (xy);
}

INVERSE(e_inverse) /* ellipsoid & spheroid */;
    if( (P->rho = hypot(xy.x, xy.y = P->rho0 - xy.y)) != 0.0 ) {
        if (P->n < 0.) {
            P->rho = -P->rho;
            xy.x = -xy.x;
            xy.y = -xy.y;
        }
        lp.phi =  P->rho / P->dd;
        if (P->ellips) {
            lp.phi = (P->c - lp.phi * lp.phi) / P->n;
            if (fabs(P->ec - fabs(lp.phi)) > TOL7) {
                if ((lp.phi = phi1_(lp.phi, P->e, P->one_es)) == HUGE_VAL)
                    I_ERROR
            } else
                lp.phi = lp.phi < 0. ? -HALFPI : HALFPI;
        } else if (fabs(lp.phi = (P->c - lp.phi * lp.phi) / P->n2) <= 1.)
            lp.phi = asin(lp.phi);
        else
            lp.phi = lp.phi < 0. ? -HALFPI : HALFPI;
        lp.lam = atan2(xy.x, xy.y) / P->n;
    } else {
        lp.lam = 0.;
        lp.phi = P->n > 0. ? HALFPI : - HALFPI;
    }
    return (lp);
}

So the define should then become:

#define PROJ_PARMS__ \
    double  ec; \
    double  n; \
    double  c; \
    double  dd; \
    double  n2; \
    double  rho0; \
    double  phi1; \
    double  phi2; \
    double  *en; \
    int     ellips;

So, change the above to:

FORWARD(e_forward); /* ellipsoid & spheroid */
    double rho;
    if ((rho = P->c - (P->ellips ? P->n * pj_qsfn(sin(lp.phi),
        P->e, P->one_es) : P->n2 * sin(lp.phi))) < 0.) F_ERROR
    rho = P->dd * sqrt(rho);
    xy.x = rho * sin( lp.lam *= P->n );
    xy.y = P->rho0 - rho * cos(lp.lam);
    return (xy);
}

INVERSE(e_inverse) /* ellipsoid & spheroid */;
    double rho;
    if( (rho = hypot(xy.x, xy.y = rho0 - xy.y)) != 0.0 ) {
        if (P->n < 0.) {
            rho = -rho;
            xy.x = -xy.x;
            xy.y = -xy.y;
        }
        lp.phi =  rho / P->dd;
        if (P->ellips) {
            lp.phi = (P->c - lp.phi * lp.phi) / P->n;
            if (fabs(P->ec - fabs(lp.phi)) > TOL7) {
                if ((lp.phi = phi1_(lp.phi, P->e, P->one_es)) == HUGE_VAL)
                    I_ERROR
            } else
                lp.phi = lp.phi < 0. ? -HALFPI : HALFPI;
        } else if (fabs(lp.phi = (P->c - lp.phi * lp.phi) / P->n2) <= 1.)
            lp.phi = asin(lp.phi);
        else
            lp.phi = lp.phi < 0. ? -HALFPI : HALFPI;
        lp.lam = atan2(xy.x, xy.y) / P->n;
    } else {
        lp.lam = 0.;
        lp.phi = P->n > 0. ? HALFPI : - HALFPI;
    }
    return (lp);
}

[rouault]

Your suggestion looks reasonable at first, but I'm not sure that the forward and inverse methods have been intended to be thread-safe by using the same PJ* structure (or if they were, the reality is that they are not). I'd suggest rather to instanciate a PJ* structure per thread. There are for example a few projections that call pj_ctx_set_errno() (healpix, tmerc, ...) in case of errors, and the wrappers pj_fwd() & pj_inv() also do it, so this cannot be thread safe. And looking closely in the above code, I see that the I_ERROR macro also calls pj_ctx_set_errno()...

So we could apply your change, but mainly for the sake of the clarity of the code, and not as a thread-safety guarantee (by the way a proper pull request would be appreciated to ease the work of us, lazy maintainers). @hobu any opinion ?

[hobu]

What's the emoji for sticking my fingers in my ears and going "la la la la"?

[QuLogic]

How's :hear_no_evil:?

[hobu]

How's :hear_no_evil:?

That's the one!

Anyway, I'm hesitant to twist up the code too much in the name of thread safety unless someone is making that a major priority and going through stuff to tighten it up and back it up with tests. As @rouault stated, this patch fixes an immediate issue, but further issues still lurk.