Opened 16 years ago

Last modified 15 years ago

#811 closed defect (fixed)

NTF / Gradians handling issue?

Reported by: warmerdam Owned by: warmerdam
Priority: high Milestone:
Component: OGR_SRS Version: unspecified
Severity: normal Keywords:
Cc: hsaggu@…

Description

Hi Frank,
We were investigating a bug report from user where some of the EPSG NTF
coordinate system definitions have a unit error. Here is the description
from the bug report.
------------------------------------------------------------------------
Many of the EPSG NTF coordinate system definitions in FME (2004 ICE or
Beta 1521) have major errors. It appears that Gradients have been
confused with degrees. For example ESPG:27562 properties shows a
latitude origin of 52 degrees when it should be 52 Gradients or 46.8
degrees. This error seems to be replicated across many other EPSG NTF
definitions in the Coord Sys gallery
(EPSG:27561,62,63,71,72,73,81,82,83,91,92,93).
------------------------------------------------------------------------
-

The problem may be in the following method of ogrspatialreference.cpp

OGRErr OGRSpatialReference::SetNormProjParm( const char * pszName,
                                            double dfValue )

{
   GetNormInfo();

   if( (dfToDegrees != 1.0 || dfFromGreenwich != 0.0)
       && IsAngularParameter(pszName) )
   {
#ifdef WKT_LONGITUDE_RELATIVE_TO_PM
       if( dfFromGreenwich != 0.0 && IsLongitudeParameter( pszName ) )
           dfValue -= dfFromGreenwich;
#endif

       // Frank should check this, but it seems that for ESPG 27561
(for example)
       // the French ones, the anglular measure was in gradients.
However, when
       // EPSGGetProjTRFInfo was called, it already put the parameter
in as degrees,
       // a conversion was done at that time. Doing the conversion
again puts us
       // back to graidents.  It seems we should do nothing here.
       // (Dale/Harbinder 20050316)

       // dfValue /= dfToDegrees;
   }
   else if( dfToMeter != 1.0 && IsLinearParameter( pszName ) )
       dfValue /= dfToMeter;

   return SetProjParm( pszName, dfValue );
}

Please check the commented out section and the line "// dfValue /=
dfToDegrees;"

If you need more information please let me know.

Regards,
Harbinder

Change History (3)

comment:1 Changed 15 years ago by mapserver@…

I encounter a problem in the reprojection of EPSG 27572 to EPSG 4326.
My reprojected layer are nearly 400 km on the left of the other layer.
I use gdal 1.3.1 and proj 4.4.9

When i manually change the original epsg entry 
(
<27572> +proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=-2.33722917
+k_0=0.99987742 +x_0=600000 +y_0=2200000 +a=6378249.2 +b=63565
15 +towgs84=-168,-60,320,0,0,0,0 +pm=paris +units=m +no_defs  <>
)
which comes with proj 4.4.9 for EPSG 27572 to:

<27572> +proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=2.33722917
+k_0=0.99987742 +x_0=600000 +y_0=2200000 +a=6378249.2 +b=6356515
+towgs84=-168,-60,320,0,0,0,0 +units=m +no_defs  <>

or

<27572> +proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=0 +k_0=0.99987742
+x_0=600000 +y_0=2200000 +a=6378249.2 +b=6356515
+towgs84=-168,-60,320,0,0,0,0 +pm=paris +units=m +no_defs  <>

the shift disappear.

I have no idea what happens behind the curtain but i assume that it has
something to do with the prime meridian. My be a sign error.

comment:2 Changed 15 years ago by warmerdam

I have reviewed this problem, and it seems like longitudinal parameters
for projections in EPSG projections are supposed to be relative to greenwich
as is the central meridian for EPSG:26591 (stored as 9 degrees east of greenwich)
even with using a different prime meridian. 

But the longitudes for NTF seem to be stored relative to the paris meridian. 

Interestingly the record for coord_op_code = 18092 (the operation code
for EPSG:27562) has the remark field "Longitude is referenced to the Paris
meridian. Superseded in 1972 by Lambert zone II (code 18082).".  To me
this seems to suggest that storing longitudes relative to the prime meridian of
the geogcs is not typical. 

I searched in vain for some indication in the projection parameters what
prime meridian they should be taken relative to, but I couldn't find any
pertinent differences between EPSG:26591 and EPSG:27562. 

For now, I don't see how to resolve this problem properly, so I am just going
to leave it open. 

I will follow up with EPSG and try to remember to append any responses here. 

comment:3 Changed 15 years ago by warmerdam

OK, I have established that the EPSG longitude parameters are actually
relative to the geogcs of the coorinate systems's meridian. 

So I have modified ogr_srs_proj4.c accordingly and added tests for prime
meridian support in gdalautotest/osr/{osr_pm.py,osr_ct_proj4.cpp}.

Note: See TracTickets for help on using tickets.