Re: [PROJ] Create a custom CRS using EPSG transform method 9666 (P6 I=J+90 seismic bin grid coordinate operation) with pyproj

[Even Rouault via PROJ]

Yannick,

While PROJ can parse any valid WKT (by the way the SCALEUNIT[] inside 
the AXIS[] are normally not allowed for a ordinal CS, but that's not 
critical here), it can't necessarily always make practical use of it. In 
that instance, the underlying maths for `P6 I=J+90 seismic bin grid 
coordinate operation` (EPSG code `9666`) is not currently implemented in 
PROJ. Quickly looking at the formulas documented in the EPSG guidance 
note 7-2, it looks at first glance to boil down to just an affine 
transformation once you've compacted all the formulas, so you could 
potentially use a PROJ pipeline using the +proj=affine operator, or 
rewrite your WKT using the formalism of the Affine Parametric 
Transformation (EPSG:9624) implemented by PROJ.

Even

Le 02/03/2024 à 15:45, Yannick Portier via PROJ a écrit :
Hello all,

I tried asking this question on SO but failed to get an answer, so 
here it is again:

I am trying to recreate the [EPSG example of georeferencing I=J+90 
local bin 
grid](https://epsg.org/transformation_6918/EPSG-example-of-georeferencing-I-J-90-local-bin-grid.html) 
using the `P6 I=J+90 seismic bin grid coordinate operation` (EPSG code 
`9666`) but the resulting CRS can't be used to transform coordinates.
`xform = Transformer.from_crs('EPSG:25832', my_crs)` fails with 
`ProjError: Input is not a transformation.`

Here is the `wkt` string I'm using to create the crs with `my_crs = 
crs.CRS.from_wkt(wkt)`.

    DERIVEDPROJCRS["seismic survey bin grid",
        BASEPROJCRS["WGS 84 / UTM zone 31N",
            BASEGEOGCRS["WGS 84",
                ENSEMBLE["World Geodetic System 1984 ensemble",
                    MEMBER["World Geodetic System 1984 (Transit)"],
                    MEMBER["World Geodetic System 1984 (G730)"],
                    MEMBER["World Geodetic System 1984 (G873)"],
                    MEMBER["World Geodetic System 1984 (G1150)"],
                    MEMBER["World Geodetic System 1984 (G1674)"],
                    MEMBER["World Geodetic System 1984 (G1762)"],
                    MEMBER["World Geodetic System 1984 (G2139)"],
                    ELLIPSOID["WGS 84",6378137,298.257223563,
                        LENGTHUNIT["metre",1]],
                    ENSEMBLEACCURACY[2.0]],
                PRIMEM["Greenwich",0,
                    ANGLEUNIT["degree",0.0174532925199433]]],
            CONVERSION["UTM zone 31N",
                METHOD["Transverse Mercator",
                    ID["EPSG",9807]],
                PARAMETER["Latitude of natural origin",0,
                    ANGLEUNIT["degree",0.0174532925199433],
                    ID["EPSG",8801]],
                PARAMETER["Longitude of natural origin",3,
                    ANGLEUNIT["degree",0.0174532925199433],
                    ID["EPSG",8802]],
                PARAMETER["Scale factor at natural origin",0.9996,
                    SCALEUNIT["unity",1],
                    ID["EPSG",8805]],
                PARAMETER["False easting",500000,
                    LENGTHUNIT["metre",1],
                    ID["EPSG",8806]],
                PARAMETER["False northing",0,
                    LENGTHUNIT["metre",1],
                    ID["EPSG",8807]]]],
        DERIVINGCONVERSION["seismic survey bin grid",
            METHOD["P6 I=J+90 seismic bin grid coordinate operation",
                ID["EPSG",9666]],
            PARAMETER["Bin grid origin I",1,
                SCALEUNIT["Bin",1],
                ID["EPSG",8733]],
            PARAMETER["Bin grid origin J",1,
                SCALEUNIT["Bin",1],
                ID["EPSG",8734]],
            PARAMETER["Bin grid origin Easting",456781,
                LENGTHUNIT["metre",1],
                ID["EPSG",8735]],
            PARAMETER["Bin grid origin Northing",5836723,
                LENGTHUNIT["metre",1],
                ID["EPSG",8736]],
            PARAMETER["Scale factor of bin grid",0.99984,
                SCALEUNIT["Unity",1],
                ID["EPSG",8737]],
            PARAMETER["Bin width on I-axis",25,
                LENGTHUNIT["metre",1],
                ID["EPSG",8738]],
            PARAMETER["Bin width on J-axis",12.5,
                LENGTHUNIT["metre",1],
                ID["EPSG",8739]],
            PARAMETER["Map grid bearing of bin grid J-axis",20,
                ANGLEUNIT["degree",0.0174532925199433],
                ID["EPSG",8740]],
            PARAMETER["Bin node increment on I-axis",1,
                SCALEUNIT["Bin",1],
                ID["EPSG",8741]],
            PARAMETER["Bin node increment on J-axis",1,
                SCALEUNIT["Bin",1],
                ID["EPSG",8742]]],
        CS[ordinal,2],
            AXIS["inline (I)",eastSouthEast,
                ORDER[1],
                SCALEUNIT["Bin",1]],
            AXIS["crossline (J)",northNorthEast,
                ORDER[2],
                SCALEUNIT["Bin",1]],
        USAGE[
            SCOPE["Navigation and medium accuracy spatial referencing."],
            AREA["Between 0°E and 6°E, northern hemisphere between 
equator and 84°N, onshore and offshore. Algeria. Andorra. Belgium. 
Benin. Burkina Faso. Denmark - North Sea. France. Germany - North Sea. 
Ghana. Luxembourg. Mali. Netherlands. Niger. Nigeria. Norway. Spain. 
Togo. United Kingdom (UK) - North Sea."],
            BBOX[0,0,84,6]]]

What's wrong ?

Many thanks,
YeO

_______________________________________________
PROJ mailing list
PROJ@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/proj

--
http://www.spatialys.com
My software is free, but my time generally not.

_______________________________________________
PROJ mailing list
PROJ@lists.osgeo.org
https://lists.osgeo.org/mailman/listinfo/proj