PROJ.4 parameters for LCC (2SP) misinterpreted

[tomyun]

A PROJ.4 syntax for Lambert Conic Conformal with 2 standard parallels seems to be incorrectly parsed as 1 standard parallel parameters.

srs = osr.SpatialReference()
srs.ImportFromProj4('+proj=lcc +lat_1=30 +lat_2=60 +lat_0=30 +lon_0=126 +x_0=0 +y_0=0 +a=6378160 +b=6356775 +units=m +no_defs')
srs.ExportToProj4()

At this point, srs shows different parameters from what it has originally imported.

'+proj=lcc +lat_1=30 +lat_0=30 +lon_0=126 +k_0=1 +x_0=0 +y_0=0 +a=6378160 +b=6356775 +units=m +no_defs '

However, if we use a bit verbose WKT syntax for the same projection,

srs.ImportFromWkt('''
PROJCS["unnamed",
    GEOGCS["Coordinate System imported from GRIB file",
        DATUM["unknown",
            SPHEROID["Spheroid imported from GRIB file",6378160,298.2539162964695]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Lambert_Conformal_Conic_2SP"],
    PARAMETER["standard_parallel_1",30],
    PARAMETER["standard_parallel_2",60],
    PARAMETER["latitude_of_origin",30],
    PARAMETER["central_meridian",126],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0]]
''')
srs.ExportToProj4()

it now matches with the intended parameters of LCC (2SP).

'+proj=lcc +lat_1=30 +lat_2=60 +lat_0=30 +lon_0=126 +x_0=0 +y_0=0 +a=6378160 +b=6356775 +units=m +no_defs '

[Even Rouault]

Linked tickets #3998, #4605 and #5191

[45136]

Another case is reported here: ​http://gis.stackexchange.com/questions/81633/qgis-not-interpreting-lambert-conformal-conic-projection-correctly, resulting in an offset of more than 10 kilometres.

[Kyle Shannon]

I've poked around with gdalsrsinfo and found that the proj.4 output doesn't roundtrip:

kyle@kyle-workstation:~/Desktop$ gdalsrsinfo -o proj4 epsg:3979
'+proj=lcc +lat_1=49 +lat_2=77 +lat_0=49 +lon_0=-95 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs '

kyle@kyle-workstation:~/Desktop$ gdalsrsinfo '+proj=lcc +lat_1=49 +lat_2=77 +lat_0=49 +lon_0=-95 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs '

PROJ.4 : '+proj=lcc +lat_1=49 +lat_0=49 +lon_0=-95 +k_0=1 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs '

OGC WKT :
PROJCS["unnamed",
    GEOGCS["GRS 1980(IUGG, 1980)",
        DATUM["unknown",
            SPHEROID["GRS80",6378137,298.257222101],
            TOWGS84[0,0,0,0,0,0,0]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Lambert_Conformal_Conic_1SP"],
    PARAMETER["latitude_of_origin",49],
    PARAMETER["central_meridian",-95],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]

As you can see, the other standard parallel gets eaten when the parser automatically sets the projection to LCC 1SP because the latitude_ of_origin and standard_parallel_1 are equal.

Poking around some websites I found:

...
Options
 
Specify the first standard parallel and second standard parallel to tailor the projection to the area to be mapped.
 
Specify the latitude origin and longitude origin to center the map projection to the area to be mapped. Specifying a non-Equatorial or non-polar origin causes an oblique projection.
 
The Lambert Conformal Conic (Single Parallel) choice is a simplified version that uses only one parallel, the central latitude. This has slightly more distortion than the full version since the cone of projection is tangent to the Earth at the one parallel.

at: ​http://www.georeference.org/doc/lambert_conformal_conic.htm

so it's hard for me to tell what exactly is 'right'.

Since we interpret EPSG:3979 as a LCC 2SP, I would think it should be parsed that way. I seem to remember having to do a similar thing in the netCDF driver as well. I looked at the code briefly and I think I can supply a patch for testing, although I am leery or SRS stuff I don't fully understand.

[Kyle Shannon]

I added a simple test to the osr_proj4 round-tripping test where 2 standard parallels were defined with lat_1==latitude_of_origin. It failed. I added a simple fix to ogr_srs_proj4.cpp as a fix. Test and osr patch attached. Any testing appreciated.

[saultdon]

I'm having troubles applying the 5191_test.diff patch to GDALs 1.10.1 src before compiling.

The osr folder doesn't exist in the tar.gz package from the gdal site but I found the rev mentioned in the patch at http://trac.osgeo.org/gdal/browser/trunk/autotest?rev=26758.

Should I build from the svn source (http://trac.osgeo.org/gdal/browser/tags/1.10.1) instead of the tar.gz?

[Kyle Shannon]

Replying to saultdon:

I'm having troubles applying the 5191_test.diff patch to GDALs 1.10.1 src before compiling.

The osr folder doesn't exist in the tar.gz package from the gdal site but I found the rev mentioned in the patch at http://trac.osgeo.org/gdal/browser/trunk/autotest?rev=26758.

Should I build from the svn source (http://trac.osgeo.org/gdal/browser/tags/1.10.1) instead of the tar.gz?

The test patch is for the autotest directory. It is for regression testing, and I wouldn't worry too much about it. If you can apply the 5191.diff, and test your datasets/test cases against it, that'd be fantastic.

[saultdon]

My datasets are lining up now on linux, but in QGIS they come in assigned as

USER:100003 -  * Generated CRS (+proj=lcc +lat_1=49 +lat_2=77 +lat_0=49 +lon_0=-95 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs)

instead of EPSG:3978 or 3979 - those are still missing from the projections list in QGIS, but not a big issue since the "USER" assigned one is correct.

ogr2ogr works though when reprojecting a dataset into

-t_srs "EPSG:3978"

.

And this issue isn't present in Windows 8.1 with QGIS 2.0.1 from the 64bit OSGeo4W installer. QGIS can find the EPSG:3978 in it's list on that OS.

[Kyle Shannon]

Replying to saultdon:

My datasets are lining up now on linux, but in QGIS they come in assigned as

USER:100003 -  * Generated CRS (+proj=lcc +lat_1=49 +lat_2=77 +lat_0=49 +lon_0=-95 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs)

instead of EPSG:3978 or 3979 - those are still missing from the projections list in QGIS, but not a big issue since the "USER" assigned one is correct.

ogr2ogr works though when reprojecting a dataset into

-t_srs "EPSG:3978"

.

And this issue isn't present in Windows 8.1 with QGIS 2.0.1 from the 64bit OSGeo4W installer. QGIS can find the EPSG:3978 in it's list on that OS.

The EPSG data was updated recently (1.10.0 or 1.10.1, not sure) and may account for the differences in recognizing the EPSG code. It is recognized on my trunk build. Are your versions the same for both OS's?

[saultdon]

Replying to kyle:

The EPSG data was updated recently (1.10.0 or 1.10.1, not sure) and may account for the differences in recognizing the EPSG code. It is recognized on my trunk build. Are your versions the same for both OS's?

On Windows, through 64bit OSGeo4W installer 1.10.1 is installed, but QGIS 2.0.1 is compiled against 1.10.0 (​http://trac.osgeo.org/osgeo4w/ticket/393).

Arch Linux shows QGIS 2.0.1 compiled against 1.10.1 which is also installed. Both are compiled from my own PKGBUILDs. ​http://pastebin.com/xRq32bVr ​http://pastebin.com/t7cpif4F

Looks like both OSs have the same GDAL 1.10.1 installed.

Happy New Year.

[Kyle Shannon]

Added a check for existence of lat_2 if lat_0==lat_1 in trunk, r26827 and a test in r26828.

[Kyle Shannon]

Closing as per last comment. Please reopen if needed.