[PROJ] Create a custom CRS using EPSG transform method 9666 (P6 I=J+90 seismic bin grid coordinate operation) with pyproj
Yannick Portier
yportiercgg at gmail.com
Sun Mar 3 11:55:33 PST 2024
Hi Even,
Indeed it's just an affine projection and I can define it with EPSG 9624 -
that's currently how I'm doing it actually.
It would have been nice to define it straight with EPSG:9666 but I guess
seismic bin grids are a bit of a niche and may not be worth the effort
implementing.
At any rate, many thanks for your feedback,
Yannick
On Sat, Mar 2, 2024 at 4:12 PM Even Rouault <even.rouault at spatialys.com>
wrote:
> 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.
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20240303/fb6b2806/attachment.htm>
More information about the PROJ
mailing list