Opened 4 years ago

Closed 3 years ago

#5191 closed defect (fixed)

PROJ.4 parameters for LCC (2SP) misinterpreted

Reported by: tomyun Owned by: Kyle Shannon
Priority: normal Milestone: 1.11.0
Component: OGR_SRS Version: 1.10.0
Severity: normal Keywords:
Cc:

Description

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 '

Attachments (2)

5191.diff (552 bytes) - added by Kyle Shannon 4 years ago.
5191_test.diff (968 bytes) - added by Kyle Shannon 4 years ago.

Download all attachments as: .zip

Change History (13)

comment:1 Changed 4 years ago by Even Rouault

Component: PythonBindingsOGR_SRS
Owner: changed from hobu to warmerdam

Linked tickets #3998, #4605 and #5191

comment:2 Changed 4 years ago by 45136

Milestone: 1.11.0
Severity: normalmajor

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.

comment:3 Changed 4 years ago by Kyle Shannon

Owner: changed from warmerdam to Kyle Shannon
Severity: majornormal
Status: newassigned

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.

comment:4 Changed 4 years ago by 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.

Changed 4 years ago by Kyle Shannon

Attachment: 5191.diff added

Changed 4 years ago by Kyle Shannon

Attachment: 5191_test.diff added

comment:5 Changed 4 years ago by 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?

comment:6 in reply to:  5 Changed 4 years ago by 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.

comment:7 Changed 4 years ago by 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.

comment:8 in reply to:  7 ; Changed 4 years ago by 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?

comment:9 in reply to:  8 Changed 4 years ago by 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.

comment:10 Changed 4 years ago by Kyle Shannon

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

comment:11 Changed 3 years ago by Kyle Shannon

Resolution: fixed
Status: assignedclosed

Closing as per last comment. Please reopen if needed.

Note: See TracTickets for help on using tickets.