ogr2ogr lcc transformation bug

[tsw]

The following reprojection will not work correctly:

ogr2ogr -s_srs "+proj=lcc +lat_1=49 +lat_2=77 +lat_0=49 +lon_0=-95 +k_0=1 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs" -t_srs "+proj=lcc +lat_1=50 +lat_2=70 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs" -f "ESRI Shapefile" reproj.shp source.shp

It will also not work if +k_0=1 is removed from the source.

It appears that the +lat_2 parameter in the source is ignored resulting in increasing visual distortion from south to north. However if you do the following:

ogr2ogr -s_srs "+proj=lcc +lat_1=49 +lat_2=77 +lat_0=49 +lon_0=-95 +k_0=1 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs" -f KML intermediate.kml source.shp

and

ogr2ogr -t_srs "+proj=lcc +lat_1=50 +lat_2=70 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs" -f "ESRI Shapefile" reproj.shp intermediate.kml

Then the reprojection will work. Attempting to use this intermediate step with another format such as GML for example does not work, leading me to suspect that there is something unique in the code base of conversion to KML that causes correct parsing of the source srs that is not applied in other cases.

The same error seems to occur in using gdal to reproject a raster going the other direction as shown above. I have not tried using an intermediate step in that case.

If needed please contact me and I can try to obtain permission to share a simplified version of the source files for testing.

[Even Rouault]

Using KML as intermediate format is equivalent to specifying -t_srs EPSG:4326 with other formats. But I don't see any difference when doing direct conversion or when doing in 2 steps with EPSG:4326 as intermediate step. See with below sample point :

$ gdaltransform  -s_srs "+proj=lcc +lat_1=49 +lat_2=77 +lat_0=49 +lon_0=-95 +k_0=1 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs" -t_srs "+proj=lcc +lat_1=50 +lat_2=70 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"
34161.6137934215 1353754.47982796
85990.8073219527 2351636.06714891 0
^C
$ gdaltransform  -s_srs "+proj=lcc +lat_1=49 +lat_2=77 +lat_0=49 +lon_0=-95 +k_0=1 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs" -t_srs EPSG:4326
34161.6137934215 1353754.47982796
-94.382708447192 61.0628801693071 0
^C
$ gdaltransform -s_srs EPSG:4326 -t_srs  "+proj=lcc +lat_1=50 +lat_2=70 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"
-94.382708447192 61.0628801693071 
85990.8073219542 2351636.06714891 0

Also note taht as far as "+proj=lcc +lat_1=49 +lat_2=77 +lat_0=49 +lon_0=-95 +k_0=1 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs" is concerned, when +lat_1 is equal to +lat_0, the code that parses the proj4 string explicitely ignores +lat_2 to use the LCC 1SP form of the projection. I'm not versed enough in projections to comment more about that, but it seems to be written voluntary in that way. And using EPSG:4326 as intermediate format will certainly not make +lat_2 be taken into account.

[warmerdam]

I'm a bit lost in this, but I will observe that I don't think it makes sense to have +lat_0, +lat_1 and +lat_2 for LCC. It appears this is being treated as LCC 1SP instead of LCC 2SP.

eg.

warmerda@gdal65[110]% gdalsrsinfo '+proj=lcc +lat_1=49 +lat_2=77 +lat_0=49 +lon_0=-95 +k_0=1 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +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["NAD83",
        DATUM["North_American_Datum_1983",
            SPHEROID["GRS 1980",6378137,298.257222101,
                AUTHORITY["EPSG","7019"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6269"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9108"]],
        AUTHORITY["EPSG","4269"]],
    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]]

So far, I see no indication of an error, just limitations in detecting and reporting odd inputs.

[tsw]

I've added a .shp file and a text file outlining the various results demonstrating the problem.

[tsw]

Thanks for your replies.

I've now attached a file with a source .shp file and the commands I've used. You will see that the results of the GLM with EPSG 4326 specified don't match the KML output (viewing in QGIS). It appears however upon further examination that my specification of a source projection does not entirely over-ride the .prj file because if it is removed I get a warning of :

ERROR 1: Latitude 1361824.629630 is invalid. Valid range is [-90,90]. This warning will not be issued any more Warning 1: Longitude 184611.241353 has been modified to fit into range [-180,180]. This warning will not be issued any more

I realize that having lat_1 and lat_0 equal doesn't match the correct specification of lcc, however, when receiving files in this format, I would like to be able to correctly convert them to a properly formatted projection.

Please review the attachment and if I have made operator errors I would appreciate knowing what I've done incorrectly.

[Even Rouault]

Linked tickets #3998, #4605 and #5191

[Kyle Shannon]

Replying to warmerdam:

I'm a bit lost in this, but I will observe that I don't think it makes sense to have +lat_0, +lat_1 and +lat_2 for LCC. It appears this is being treated as LCC 1SP instead of LCC 2SP.

eg.

warmerda@gdal65[110]% gdalsrsinfo '+proj=lcc +lat_1=49 +lat_2=77 +lat_0=49 +lon_0=-95 +k_0=1 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +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["NAD83",
        DATUM["North_American_Datum_1983",
            SPHEROID["GRS 1980",6378137,298.257222101,
                AUTHORITY["EPSG","7019"]],
            TOWGS84[0,0,0,0,0,0,0],
            AUTHORITY["EPSG","6269"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9108"]],
        AUTHORITY["EPSG","4269"]],
    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]]

So far, I see no indication of an error, just limitations in detecting and reporting odd inputs.

Frank, I've been browsing some documentation, mostly Snyder USGS publications, and I haven't seen a restriction stating that if lat_0 == lat_1 then it's LCC 1SP. It seems the reverse is true, if LCC 1SP then lat_0 == lat_1, and it's a 'tangential' form of LCC

The parser checks for lat_0 == lat_1 and if that condition is true, it's a 1SP, regardless if lat_2 is specified. The patch(es) on #5191 handle this. Gerald's note on ​http://www.remotesensing.org/geotiff/proj_list/lambert_conic_conformal_1sp.html doesn't eliminate this option either, as far as I understand.

[Even Rouault]

I'd say the fix for #5191 fixed that one.