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

[Even Rouault]

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 at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/proj

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