[Proj] mercator projection to spherical ellipsoid from WGS84

[Frank Warmerdam]

Scott Ellington wrote:
> Hi Frank,
> 
> Thanks for responding.  Unfortunately, the issue manifests itself in
> mapserver so I don't think I can chain two commands together.
> 
> As for revisiting this issue, I think this would be helpful.  At the very
> least, I'd like the ability to avoid the datum transform as it used to be
> possible to do.  I can file a bug to keep this on the radar if you'd like.

Scott,

I see what you mean.  OK, file a ticket against PROJ.4.

I'm interested in suggestions from the community on what pj_datum_transform()
ought to do when the source and destination coordinate systems have different
ellipsoids, but one, the other, or both lack a definition of how to get
to WGS84.

Currently the code looks like this, which I hope is reasonably understandable
based on function names and comments.

int pj_datum_transform( PJ *srcdefn, PJ *dstdefn,
                         long point_count, int point_offset,
                         double *x, double *y, double *z )

{
     double      src_a, src_es, dst_a, dst_es;
     int         z_is_temp = FALSE;

     pj_errno = 0;

/* -------------------------------------------------------------------- */
/*      Short cut if the datums are identical.                          */
/* -------------------------------------------------------------------- */
     if( pj_compare_datums( srcdefn, dstdefn ) )
         return 0;

     src_a = srcdefn->a;
     src_es = srcdefn->es;

     dst_a = dstdefn->a;
     dst_es = dstdefn->es;

/* -------------------------------------------------------------------- */
/*      Create a temporary Z array if one is not provided.              */
/* -------------------------------------------------------------------- */
     if( z == NULL )
     {
         int	bytes = sizeof(double) * point_count * point_offset;
         z = (double *) pj_malloc(bytes);
         memset( z, 0, bytes );
         z_is_temp = TRUE;
     }

#define CHECK_RETURN {if( pj_errno != 0 && (pj_errno > 0 || 
transient_error[-pj_errno] == 0) ) { if( z_is_temp ) pj_dalloc(z); return 
pj_errno; }}

/* -------------------------------------------------------------------- */
/*	If this datum requires grid shifts, then apply it to geodetic   */
/*      coordinates.                                                    */
/* -------------------------------------------------------------------- */
     if( srcdefn->datum_type == PJD_GRIDSHIFT )
     {
         pj_apply_gridshift( pj_param(srcdefn->params,"snadgrids").s, 0,
                             point_count, point_offset, x, y, z );
         CHECK_RETURN;

         src_a = SRS_WGS84_SEMIMAJOR;
         src_es = SRS_WGS84_ESQUARED;
     }

     if( dstdefn->datum_type == PJD_GRIDSHIFT )
     {
         dst_a = SRS_WGS84_SEMIMAJOR;
         dst_es = SRS_WGS84_ESQUARED;
     }

/* ==================================================================== */
/*      Do we need to go through geocentric coordinates?                */
/* ==================================================================== */
     if( src_es != dst_es || src_a != dst_a
         || srcdefn->datum_type == PJD_3PARAM
         || srcdefn->datum_type == PJD_7PARAM
         || dstdefn->datum_type == PJD_3PARAM
         || dstdefn->datum_type == PJD_7PARAM)
     {
/* -------------------------------------------------------------------- */
/*      Convert to geocentric coordinates.                              */
/* -------------------------------------------------------------------- */
         pj_geodetic_to_geocentric( src_a, src_es,
                                    point_count, point_offset, x, y, z );
         CHECK_RETURN;

/* -------------------------------------------------------------------- */
/*      Convert between datums.                                         */
/* -------------------------------------------------------------------- */
         if( srcdefn->datum_type == PJD_3PARAM
             || srcdefn->datum_type == PJD_7PARAM )
         {
             pj_geocentric_to_wgs84( srcdefn, point_count, point_offset,x,y,z);
             CHECK_RETURN;
         }

         if( dstdefn->datum_type == PJD_3PARAM
             || dstdefn->datum_type == PJD_7PARAM )
         {
             pj_geocentric_from_wgs84( dstdefn, point_count,point_offset,x,y,z);
             CHECK_RETURN;
         }

/* -------------------------------------------------------------------- */
/*      Convert back to geodetic coordinates.                           */
/* -------------------------------------------------------------------- */
         pj_geocentric_to_geodetic( dst_a, dst_es,
                                    point_count, point_offset, x, y, z );
         CHECK_RETURN;
     }

/* -------------------------------------------------------------------- */
/*      Apply grid shift to destination if required.                    */
/* -------------------------------------------------------------------- */
     if( dstdefn->datum_type == PJD_GRIDSHIFT )
     {
         pj_apply_gridshift( pj_param(dstdefn->params,"snadgrids").s, 1,
                             point_count, point_offset, x, y, z );
         CHECK_RETURN;
     }

     if( z_is_temp )
         pj_dalloc( z );

     return 0;
}

Best regards,
-- 
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up   | Frank Warmerdam, warmerdam at pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush    | President OSGeo, http://osgeo.org