Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Equal Area Cylindrical not using k_0 scale factor in WKT form #1881

Closed
snowman2 opened this issue Jan 25, 2020 · 6 comments · Fixed by #1884
Closed

Equal Area Cylindrical not using k_0 scale factor in WKT form #1881

snowman2 opened this issue Jan 25, 2020 · 6 comments · Fixed by #1884
Labels

Comments

@snowman2
Copy link
Contributor

@snowman2 snowman2 commented Jan 25, 2020

Example of problem

This is expected:

$ ./projinfo '+proj=cea +type=crs'
PROJ.4 string:
+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +type=crs

WKT2:2019 string:
PROJCRS["unknown",
    BASEGEOGCRS["unknown",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8901]]],
    CONVERSION["unknown",
        METHOD["Lambert Cylindrical Equal Area",
            ID["EPSG",9835]],
        PARAMETER["Latitude of 1st standard parallel",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8823]],
        PARAMETER["Longitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]

This is missing the scaling factor:

$ ./projinfo '+proj=cea +k_0=2 +type=crs'
PROJ.4 string:
+proj=cea +k_0=2 +type=crs

WKT2:2019 string:
PROJCRS["unknown",
    BASEGEOGCRS["unknown",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8901]]],
    CONVERSION["unknown",
        METHOD["Lambert Cylindrical Equal Area",
            ID["EPSG",9835]],
        PARAMETER["Latitude of 1st standard parallel",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8823]],
        PARAMETER["Longitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]

Problem description

This defines a scaling factor. https://proj.org/operations/projections/cea.html
However, it is not being used.

Expected Output

I would expect it to have a scaling factor similar to lcc or merc:

$ ./projinfo '+proj=lcc +lat_1=2 +k_0=2 +type=crs'
PROJ.4 string:
+proj=lcc +lat_0=0 +lon_0=0 +lat_1=2 +lat_2=0 +x_0=0 +y_0=0 +k_0=2 +datum=WGS84 +units=m +no_defs +type=crs

WKT2:2019 string:
PROJCRS["unknown",
    BASEGEOGCRS["unknown",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8901]]],
    CONVERSION["unknown",
        METHOD["Lambert Conic Conformal (2SP Michigan)",
            ID["EPSG",1051]],
        PARAMETER["Latitude of false origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8821]],
        PARAMETER["Longitude of false origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8822]],
        PARAMETER["Latitude of 1st standard parallel",2,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8823]],
        PARAMETER["Latitude of 2nd standard parallel",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8824]],
        PARAMETER["Easting at false origin",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8826]],
        PARAMETER["Northing at false origin",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8827]],
        PARAMETER["Ellipsoid scaling factor",2,
            SCALEUNIT["unity",1],
            ID["EPSG",1038]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]

Environment Information

$ ./projinfo
Rel. 7.0.0, March 1st, 2020
  • Ubuntu 18.04.3 LTS

Installation method

  • From source master
@snowman2 snowman2 added the bug label Jan 25, 2020
@rouault
Copy link
Member

@rouault rouault commented Jan 25, 2020

Hum this is interesting... This was ported from GDAL 2.x OSRSetCEA() to PROJ.4 method, which has no scale factor. I presume GDAL implementation comes from EPSG 'Lambert Cylindrical Equal Area' method itself (https://www.epsg-registry.org/export.htm?gml=urn:ogc:def:method:EPSG::9835) that doesn't define a scale factor. Hence this is lost when mapping to WKT.
The clean way of solving this would probably when detecting +proj=cea and +k/k_0 to use a custom WKT representation based on the EPSG one but with an extra parameter...

@snowman2
Copy link
Contributor Author

@snowman2 snowman2 commented Jan 25, 2020

I haven't dug into the PROJ code to check, but does it use the scale factor parameter?

@rouault
Copy link
Member

@rouault rouault commented Jan 25, 2020

yes, the cea.cpp file uses P->k_0

@snowman2
Copy link
Contributor Author

@snowman2 snowman2 commented Jan 25, 2020

I added the extra parameter via WKT, but then it is ignored in the PROJ string:

./projinfo 'CONVERSION["unknown",METHOD["Lambert Cylindrical Equal Area"],PARAMETER["Latitude of 1st standard parallel",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8823]],PARAMETER["Longitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["False easting",0,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]],PARAMETER["Scale factor at natural origin",1,SCALEUNIT["unity",1],ID["EPSG",8805]]]'
PROJ string:
+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0

WKT2:2019 string:
CONVERSION["unknown",
    METHOD["Lambert Cylindrical Equal Area"],
    PARAMETER["Latitude of 1st standard parallel",0,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8823]],
    PARAMETER["Longitude of natural origin",0,
        ANGLEUNIT["degree",0.0174532925199433],
        ID["EPSG",8802]],
    PARAMETER["False easting",0,
        LENGTHUNIT["metre",1],
        ID["EPSG",8806]],
    PARAMETER["False northing",0,
        LENGTHUNIT["metre",1],
        ID["EPSG",8807]],
    PARAMETER["Scale factor at natural origin",1,
        SCALEUNIT["unity",1],
        ID["EPSG",8805]]]

@rouault
Copy link
Member

@rouault rouault commented Jan 25, 2020

except that it won't receive it now, unless you use a non-CRS PROJ pipeline string

I added the extra parameter via WKT, but then it is ignored in the PROJ string

'Expected' as there's no mapping between 'Scale factor at natural origin' and k_0 for this METHOD

@snowman2
Copy link
Contributor Author

@snowman2 snowman2 commented Jan 25, 2020

Okay, so this won't work via WKT until the fix is applied. Makes sense.

rouault added a commit to rouault/PROJ that referenced this issue Jan 25, 2020
Fixes OSGeo#1881

Digging into the implementation of proj=cea, it appears that
k_0 and lat_ts are intended to be exclusive ways of specifying the
same concept. EPSG only models the variant using lat_s.
So if k_0 is found and lat_ts is absent, compute the equivalent
value of lat_ts from k_0.

Note: k_0 should normally be in the [0,1] range. In case creative users
would use something outside, we raise an exception, even if the cea
implementation could potentially deal with any k_0 value. Hopefully
this is a (reasonable) limitation that will address nominal use cases.
rouault added a commit to rouault/PROJ that referenced this issue Feb 6, 2020
Fixes OSGeo#1881

Digging into the implementation of proj=cea, it appears that
k_0 and lat_ts are intended to be exclusive ways of specifying the
same concept. EPSG only models the variant using lat_s.
So if k_0 is found and lat_ts is absent, compute the equivalent
value of lat_ts from k_0.

Note: k_0 should normally be in the [0,1] range. In case creative users
would use something outside, we raise an exception, even if the cea
implementation could potentially deal with any k_0 value. Hopefully
this is a (reasonable) limitation that will address nominal use cases.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

Successfully merging a pull request may close this issue.

2 participants