Bug report #719

Problem with NTF projection in ECW files

Added by jrepetto - over 13 years ago. Updated about 11 years ago.

Status:Closed
Priority:Low
Assignee:nobody -
Category:Projection Support
Affected QGIS version: Regression?:No
Operating System:Gentoo Easy fix?:No
Pull Request or Patch supplied: Resolution:fixed
Crashes QGIS or corrupts data: Copied to github as #:10778

Description

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",
SPHEROIDCLA80IGN,
PRIMEMGreenwich,
UNITdegree,
PROJECTIONLambert_Conformal_Conic_2SP,
PARAMETERstandard_parallel_1,
PARAMETERstandard_parallel_2,
PARAMETERlatitude_of_origin,
PARAMETERcentral_meridian,
PARAMETERfalse_easting,
PARAMETERfalse_northing,
UNITMeter
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.

History

#1 Updated by jrepetto - over 13 years ago

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.

#2 Updated by jrepetto - almost 13 years ago

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.

#3 Updated by jrepetto - almost 13 years ago

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;

#4 Updated by jrepetto - almost 13 years ago

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;

#5 Updated by jrepetto - almost 13 years ago

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.

#6 Updated by Tim Sutton almost 13 years ago

  • Status changed from Open to Closed
  • Resolution set to fixed

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

#7 Updated by Anonymous about 11 years ago

Milestone Version 0.9.1 deleted

Also available in: Atom PDF