Problem with NTF projection in ECW files

[jrepetto]

I want to use Qgis to display ECW files (french maps), such as : <​http://www.aix-mrs.iufm.fr/formations/filieres/hge/gd/gdticehg/crigepaca/donneesmartigues/scan25000martigues.ecw>

The projection recorded in the ECW file header is LM2FRANC, it corresponds to NTF Lambert II (QGIS SRSID: 1605, PostGIS SRID: 27572).

When I open the file in Qgis 0.8, an error message is displayed on the console : QgsSpatialRefSys::createFromProj4 error proj string supplied has no +ellps argument

and if I check the properties of the layer, the SRS is "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs". I have to manually change the projection.

gdalinfo gives the right projection :

$ gdalinfo scan25000martigues.ecw Driver: ECW/ERMapper Compressed Wavelets Size is 1353, 1153 Coordinate System is: PROJCS["unnamed",

GEOGCS["N.T.F.",

DATUM["NTF",

SPHEROID["CLA80IGN",6378249.2,293.4660213]],

PRIMEM["Greenwich",0], UNIT["degree",0.0174532925199433]],

PROJECTIONLambert_Conformal_Conic_2SP, PARAMETER["standard_parallel_1",45.89891888888725], PARAMETER["standard_parallel_2",47.69601444444274], PARAMETER["latitude_of_origin",46.79999999999734], PARAMETER["central_meridian",2.337229169999754], PARAMETER["false_easting",600000], PARAMETER["false_northing",2200000], UNIT["Meter",1]]

Origin = (818273.000000000000000,1827901.000000000000000) Pixel Size = (2.500000000000000,-2.500000000000000) Corner Coordinates: Upper Left ( 818273.000, 1827901.000) ( 5d 1'43.15"E, 43d25'10.72"N) Lower Left ( 818273.000, 1825018.500) ( 5d 1'38.78"E, 43d23'37.52"N) Upper Right ( 821655.500, 1827901.000) ( 5d 4'13.18"E, 43d25'6.94"N) Lower Right ( 821655.500, 1825018.500) ( 5d 4'8.74"E, 43d23'33.75"N) Center ( 819964.250, 1826459.750) ( 5d 2'55.96"E, 43d24'22.24"N) Band 1 Block=1353x1 Type=Byte, ColorInterp=Red Overviews: arbitrary Band 2 Block=1353x1 Type=Byte, ColorInterp=Green Overviews: arbitrary Band 3 Block=1353x1 Type=Byte, ColorInterp=Blue Overviews: arbitrary

I am using Qgis 0.8, gdal-1.4.1, libecwj2-3.3 and proj-4.5.0 on a Gentoo Linux system.

[jrepetto]

The missing string in the projection parameters is : +ellps=clark80 But the projection parameters in Qgis 0.8.0 have another problem (see ticket #598). Either lon_0 should be non null, either the primer meridian should be set to Paris, but not both.

[jrepetto]

As the milestone for this bug resolution is postponed to the next release whenever there is a new GIS release, I have decided to solve it myself, because it is very boring to have to set manually the projection parameters each time you open a file.

The error message is generated by the function QgsSpatialRefSys::createFromProj4 in the file qgsspatialrefsys.cpp :

QRegExp myEllipseRegExp( "\\+ellps=\\S+" );
  myStart= 0;
  myLength=0;
  myStart = myEllipseRegExp.search(theProj4String, myStart);
  if (myStart==-1)
  {
    QgsLogger::warning("QgsSpatialRefSys::createFromProj4 error proj string supplied has no +ellps argument");

    return mIsValidFlag;
  }
  else
  {
    myLength = myEllipseRegExp.matchedLength();
  }

According to the PROJ.4 documentation (available at ​ftp://ftp.remotesensing.org/proj/OF90-284.pdf), the +ellps parameter is not mandatory (see page 9, paragraph "Specifying the Earth’s Figure").

The "+ellps" parameter is only a convenient method of specifying standard ellisoidal constants.

The standard way is to use two constants (only one for a sphere) : The first and required value +a=a where a is the semimajor axis of the ellipse or equitorial radius. The second parameter can be in any one of the following standard forms:

• semiminor axis or polar radius +b=b,
  
• flattening +f=f ,
  
• reciprocal flattening, +rf=1/f ,
  
• eccentricity +e=e, or
  
• eccentricity squared +es=e2 .

In the case reported upper, the PROJ4 string is :

+proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=2.337229166666664 +k_0=0.99987742 +x_0=600000 +y_0=2200000 +a=6378249.2 +b=6356515.000000472 +units=m +no_defs

It uses the standard form to specify the ellipsoid, so QGIS should accept it.

[jrepetto]

Simple patch proposal : remove the obligation for the +ellps parameter.

diff -ur qgis_0.9.0.orig/src/core/qgsspatialrefsys.cpp qgis_0.9.0/src/core/qgsspatialrefsys.cpp
--- qgis_0.9.0.orig/src/core/qgsspatialrefsys.cpp       2007-09-16 04:45:42.000000000 +0200
+++ qgis_0.9.0/src/core/qgsspatialrefsys.cpp    2007-11-29 12:48:01.000000000 +0100
@@ -514,17 +514,11 @@
   myStart= 0;
   myLength=0;
   myStart = myEllipseRegExp.search(theProj4String, myStart);
-  if (myStart==-1)
-  {
-    QgsLogger::warning("QgsSpatialRefSys::createFromProj4 error proj string supplied has no +ellps argument");
-
-    return mIsValidFlag;
-  }
-  else
+  if (myStart!=-1)
   {
     myLength = myEllipseRegExp.matchedLength();
+    mEllipsoidAcronym = theProj4String.mid(myStart+ELLPS_PREFIX_LEN,myLength-ELLPS_PREFIX_LEN);
   }
-  mEllipsoidAcronym = theProj4String.mid(myStart+ELLPS_PREFIX_LEN,myLength-ELLPS_PREFIX_LEN);
   //mproj4string must be set here for the rest of this method to behave in a meaningful way...
   mProj4String = theProj4String;

[jrepetto]

Second patch proposal : checks that the proj4 string contain either the +ellps parameter, either the +a parameter :

diff -ur qgis_0.9.0.orig/src/core/qgsspatialrefsys.cpp qgis_0.9.0/src/core/qgsspatialrefsys.cpp
--- qgis_0.9.0.orig/src/core/qgsspatialrefsys.cpp       2007-09-16 04:45:42.000000000 +0200
+++ qgis_0.9.0/src/core/qgsspatialrefsys.cpp    2007-11-29 13:03:28.000000000 +0100
@@ -488,10 +488,13 @@
 {

   //
-  // Example:
+  // Examples:
   // +proj=tmerc +lat_0=0 +lon_0=-62 +k=0.999500 +x_0=400000 +y_0=0
   // +ellps=clrk80 +towgs84=-255,-15,71,0,0,0,0 +units=m +no_defs
   //
+  // +proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=2.337229166666664 +k_0=0.99987742
+  // +x_0=600000 +y_0=2200000 +a=6378249.2 +b=6356515.000000472 +units=m +no_defs
+  //
   mIsValidFlag=false;

   QRegExp myProjRegExp( "\\+proj=\\S+" );
@@ -514,17 +517,22 @@
   myStart= 0;
   myLength=0;
   myStart = myEllipseRegExp.search(theProj4String, myStart);
-  if (myStart==-1)
+  if (myStart!=-1)
   {
-    QgsLogger::warning("QgsSpatialRefSys::createFromProj4 error proj string supplied has no +ellps argument");
-
-    return mIsValidFlag;
+    myLength = myEllipseRegExp.matchedLength();
+    mEllipsoidAcronym = theProj4String.mid(myStart+ELLPS_PREFIX_LEN,myLength-ELLPS_PREFIX_LEN);
   }
-  else
+
+  QRegExp myAxisRegExp( "\\+a=\\S+" );
+  myStart= 0;
+  myLength=0;
+  myStart = myAxisRegExp.search(theProj4String, myStart);
+  if (myStart==-1 && mEllipsoidAcronym.isNull())
   {
-    myLength = myEllipseRegExp.matchedLength();
+    QgsLogger::warning("QgsSpatialRefSys::createFromProj4 error proj string supplied has no +ellps or +a argument");
+
+    return mIsValidFlag;
   }
-  mEllipsoidAcronym = theProj4String.mid(myStart+ELLPS_PREFIX_LEN,myLength-ELLPS_PREFIX_LEN);
   //mproj4string must be set here for the rest of this method to behave in a meaningful way...
   mProj4String = theProj4String;

[jrepetto]

Note : With one of the above patches applied, and when there is no +ellps argument, new warnings appear :

Warning: QgsSpatialRefSys::getRecord failed :  select * from tbl_srs where parameters='+proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=2.337229166666664 +k_0=0.99987742 +x_0=600000 +y_0=2200000 +a=6378249.2 +b=6356515.000000472 +units=m +no_defs'
Warning: QgsSpatialRefSys::findMatchingProj will only work if prj acr ellipsoid acr and proj4string are set!...
Warning: QgsSpatialRefSys::getRecord failed :  select * from tbl_srs where parameters='+proj=lcc +lat_1=46.8 +lat_0=46.8 +lon_0=2.337229166666664 +k_0=0.99987742 +x_0=600000 +y_0=2200000 +a=6378249.2 +b=6356515.000000472 +units=m +no_defs'
Warning: QgsSpatialRefSys::findMatchingProj will only work if prj acr ellipsoid acr and proj4string are set!...

I don't think it is a problem, the map is correctly displayed and referenced.

[timlinux]

Hi

I have applied your patch to SVN trunk as r

Please note for future patch submissions please attach them to the ticket as bug7689fix.diff as described in ​2.6. Submitting Patches.

Many thanks for your contribution! I am adding you to our bug triage page for fame & glory :-)

Regards

Tim

[(none)]

Milestone Version 0.9.1 deleted