Support for CRS described as PROJ pipeline

[jerstlouis]

Feature description

Currently it is not possible to use gdalwarp (and possibly also ogr2ogr and QGIS) to transform to/from a CRS described by a PROJ string using a pipeline.

For example,

gdalwarp  ~/Desktop/gebco.tiff test.tiff
   -s_srs "EPSG:4326"
   -t_srs "+proj=pipeline
      +step +proj=geoc
      +step +proj=isea +R=6371007.18091875 +x_0=19186144.8709401525557041 +y_0=-3323137.7717845365405083"

results in:

ERROR 1: PROJ: proj_crs_get_coordinate_system: Object is not a SingleCRS
ERROR 1: PROJ: proj_crs_get_coordinate_system: Object is not a SingleCRS
Using band 4 of destination image as alpha.
Using band 4 of source image as alpha.
Warning 1: The definition of geographic CRS EPSG:4326 got from GeoTIFF keys is not the same as the one from the EPSG registry, which may cause issues during reprojection operations. Set GTIFF_SRS_SOURCE configuration option to EPSG to use official parameters (overriding the ones from GeoTIFF keys), or to GEOKEYS to use custom values from GeoTIFF keys and drop the EPSG code.
ERROR 1: PROJ: proj_crs_get_coordinate_system: Object is not a SingleCRS
ERROR 1: PROJ: proj_crs_get_coordinate_system: Object is not a SingleCRS
ERROR 1: PROJ: proj_crs_get_datum: Object is not a SingleCRS
ERROR 1: PROJ: proj_crs_get_datum_ensemble: Object is not a SingleCRS
ERROR 1: PROJ: proj_crs_get_coordinate_system: Object is not a SingleCRS
ERROR 1: PROJ: proj_create_operations: target_crs is not a CRS or a CoordinateMetadata
ERROR 6: Cannot find coordinate operations from `EPSG:4326' to `{
  "$schema": "https://proj.org/schemas/v0.7/projjson.schema.json",
  "type": "Conversion",
  "name": "PROJ-based coordinate operation",
  "method": {
    "name": "PROJ-based operation method: +proj=pipeline +step +proj=geoc +step +proj=isea +R=6371007.18091875 +x_0=19186144.8709401525557041 +y_0=-3323137.7717845365405083 +type=crs"
  },
  "parameters": []
}

This is preventing the use of an ISEA CRS which I believe currently can only be accurately described using this pipeline, since +proj=isea does not currently automatically apply the expected geodetic -> geocentric conversion (separate issue to clarify whether that should always be done automatically).

Another use case is for the 5x6 ISEA space, where we want to perform an additional +proj=affine step.

Additional context

See also OSGeo/PROJ#4211

[rouault]

This is just a usage issue. You cannot use a pipeline as a CRS. A pipeline expresses a coordinate operation which is not a CRS. You need to use the -ct switch to specify a custom pipeline. Something like:

gdalwarp  ~/Desktop/gebco.tiff test.tiff
   -s_srs "EPSG:4326"
   -t_srs "+proj=isea +R=6371007.18091875 +x_0=19186144.8709401525557041 +y_0=-3323137.7717845365405083 +type=crs"
   -ct "+proj=pipeline +step +proj=axisswap +order=2,1 +step +proj=unitconvert +xy_in=deg +xy_out=rad +step +proj=geoc
      +step +proj=isea +R=6371007.18091875 +x_0=19186144.8709401525557041 +y_0=-3323137.7717845365405083"

[jerstlouis]

@rouault Thanks, that command does generate a test.tiff that looks right, so it might work as I intended (the difference is hard to tell from without the +proj=geoc step).

A pipeline expresses a coordinate operation which is not a CRS.

This is confusing to me, because in ISO19111, CRSes can include all types of transformation.
So I would have expected that you can define a CRS that includes any such transformation with a PROJ string.
(but then units do get confusing).

And these OGC:1534 and OGC:153456 CRSes which we want to define require this extra +proj=geoc step before the ISEA, and +proj=affine (in the case of the 5x6) steps in their definitions.

If I understand correctly, the -t_srs specifies what the target SRS will be, but when using -ct the operations to do the transformation are defined strictly by -ct (this is why the +proj=isea is in both places).

The problem with that is that currently, the +proj=isea does not automatically perform the +proj=geoc step (that might actually be a bug).

Then if accessing the test.tiff produced with that pipeline, if the CRS is defined without those extra steps, the inverse projection will result in the wrong latitude and longitude since it expects the regular +proj=isea that did not apply the +proj=geoc step.

So that +proj=geoc step really needs to be considered part of the CRS definition.

I'm really confused how such scenarios are intended to work.
Perhaps this is a limitation of PROJ strings, that PROJ JSON / CRS JSON will handle better?

[rouault]

Perhaps this is a limitation of PROJ strings, that PROJ JSON / CRS JSON will handle better?

Perhaps a solution would be to build the projected CRS on top of a geodetic CRS that uses geocentric latitudes, like done for some IAU CRS, like

PROJCRS["Mars (2015) / Ocentric / Equirectangular, clon = 0",
    BASEGEODCRS["Mars (2015) / Ocentric",
        DATUM["Mars (2015)",
            ELLIPSOID["Mars (2015)",3396190,169.894447223612,
                LENGTHUNIT["metre",1]],
            ANCHOR["Viking 1 lander: 47.95137 W"]],
        PRIMEM["Reference Meridian",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["IAU",49902,2015]],
    CONVERSION["Equirectangular, clon = 0",
        METHOD["Equidistant Cylindrical",
            ID["EPSG",1028]],
        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]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    ID["IAU",49912,2015]]

which relies on

GEODCRS["Mars (2015) / Ocentric",
    DATUM["Mars (2015)",
        ELLIPSOID["Mars (2015)",3396190,169.894447223612,
            LENGTHUNIT["metre",1]],
        ANCHOR["Viking 1 lander: 47.95137 W"]],
    PRIMEM["Reference Meridian",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[spherical,2],
        AXIS["planetocentric latitude (U)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["planetocentric longitude (V)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["IAU",49902,2015],
    REMARK["Source of IAU Coordinate systems: https://doi.org/10.1007/s10569-017-9805-5"]]

but I'm not sure that would add the +proj=geoc step the way you'd like. You're on the edge of what the modelling allows to do