[Gdal-dev] Radians as units in OGRSpatialReference and OGRProj4CT

[David E. Knight]

Hello, I have built a GIS sytem in which I use radians as my unit of measurement. To serve that purpose I created a OGRSpatialReference to represent a radian based projection.  Here is the code that I use...
////begin code snippet////
 
char* strWkt="GEOGCS[\"MY RADIAN WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",7030]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY[\"EPSG\",6326]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",8901]],UNIT[\"radian\",1.0],AXIS[\"Lat\",NORTH],AXIS[\"Long\",EAST],AUTHORITY[\"EPSG\",4326]]";
  
OGRSpatialReference myProjection;
myProjection.importFromWkt(&strWkt);
////end code snippet////
Notice that my Unit node looks like this ...
UNIT[\"radian\",1.0]
My question: is this a correct use of the wkt string?  
If the answer is no than please ignore the rest of this posting.  If the answer is yes than there seems to be a bug in the OGRProj4CT class.  OGRProj4CT is the class which handles conversions between projections.  In cases where one of the projections involved in the conversion has a radians based unit the transform method returns faulty results.  Here are the relevant bits of code...

////begin code snippet////
int OGRProj4CT::Initialize( OGRSpatialReference * poSourceIn,OGRSpatialReference * poTargetIn )
{
 poSRSSource = poSourceIn->Clone();
 poSRSTarget = poTargetIn->Clone();
 bSourceLatLong = poSRSSource->IsGeographic();
 bTargetLatLong = poSRSTarget->IsGeographic();
////end code snippet////
////begin code snippet////
int OGRProj4CT::Transform( int nCount, double *x, double *y, double *z )
{
    int   err, i;

    if( bSourceLatLong )
    {
        for( i = 0; i < nCount; i++ )
        {
            x[i] *= DEG_TO_RAD;
            y[i] *= DEG_TO_RAD;
        }
    }
//transform code snipped out//////
if( bTargetLatLong )
    {
        for( i = 0; i < nCount; i++ )
        {
            x[i] *= RAD_TO_DEG;
            y[i] *= RAD_TO_DEG;
        }
    }

////end code snippet////
As you can see the flags bSourceLatLong , bTargetLatLong  turn on conversion between radians and degrees.  The problem is that no check is made to ascertain that the unit really is degrees.  In my local codebase I have made and used the following change successfully.

////begin code snippet////
int OGRProj4CT::Initialize( OGRSpatialReference * poSourceIn,OGRSpatialReference * poTargetIn )
{
 poSRSSource = poSourceIn->Clone();
 poSRSTarget = poTargetIn->Clone();
 bSourceLatLong=FALSE;
 bTargetLatLong=FALSE;
 if(poSRSSource->IsGeographic())
 {
  const char * unitType=poSRSSource->GetAttrValue("GEOGCS|UNIT",0);
  bSourceLatLong=(strstr(unitType,"degree")!=NULL);
 }
 if(poSRSTarget->IsGeographic())
 {
  const char *unitType=poSRSTarget->GetAttrValue("GEOGCS|UNIT",0);
  bTargetLatLong=(strstr(unitType,"degree")!=NULL);
 }

////end code snippet////

My apologies if this topic has been discussed before.  The Yahoo mailing list archives are currently inaccessible. 
 
Thanks,
David E. Knight 
david1<*at*>davideknight<*dot*>com
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/gdal-dev/attachments/20021215/d455b356/attachment.html