convertTrajectory enhancements for scale and interpolation

CC @kiranallamsety @joshtownsend - - placeholder (I cannot seem to tag Phil or Roy)

Following changes are proposed to API convertTrajectory v3, which may result in a v4 change to request and response, with objectives:

1) Return point scale factor and grid convergence at each station (TBD: depth correction factor).
2) Add a new "method": "GNL" as alternative to "AzimuthalEquidistant" and "LMP", which does not apply scaling.
3) Add an input "MD_i" to the request body on which to interpolate.

Note: bug #52 (closed) should be fixed first (or at the same time as addressing these enhancements). Related issue #26 (closed) with this service can be closed when item 1 is completed.

Ad 1)

Current response has a property "stations" per api specification. Each "station" contains for example the azimuthTN and azimuthGN. The objective of this item is to add an additional output field for "point scale factor" (psf) and grid convergence at each station.

Note: There are various ways in which psf can be computed, for example we could call the Apache SIS engine at all computed trajectory (N,E) locations in the projected CRS. That would get the correctly computed actual psf. However, that approach would require an additional specific call to SIS. To keep this simple and generic, the math/method that is already implemented is used instead, and the scale factor and convergence are computed from the response as described here.

requires a projected CRS for the calculated trajectory (which is there)
For method "azimuthalEquidistant" and "LMP" the psf can be computed as described and simply added to the output as extra property.
- For the new method "GNL" (grid north local) the psf is added by the same trick. i.e., "GNL" is implemented by calling "azimuthalEquidistant" and then "unscaling" the output, see Ad 2) below.

Example desired response (note the additional "scalefactor" and "convergence" properties):

    "stations": [
        {
            "md": 0.0,
            "inclination": 10.0,
            "azimuthTN": 100.51354989131318,
            "azimuthGN": 100.0,
            "dxTN": 0.0,
            "dyTN": 0.0,
            "scalefactor": 1.000241,
            "convergence": 0.51355,
            "point": {
                "x": 1999999.99999999,
                "y": 9999999.99999969,
                "z": 0.0
            },
            "wgs84Longitude": -85.88980921169528,
            "wgs84Latitude": 27.553258329196616,
            "dls": 0.0,
            "original": true,
            "dz": 0.0
        },
...

scalefactor rounded to 6 decimal places (1 mm per km).
convergence rounded to 5 decimal places (0.2mm per km), computed as GC=TN-GN, and wrapped to interval (-180,+180) by "if convergence<-180 then convergence+=360; if convergence>180 then convergence-=360".

TBD if depth correction factor should be returned. This is a simple correction factor computed as depth/radius (see doc). For now the decision is that this is not needed (because it is easily calculated if needed).

Update 2023-05-10 BK: A mistake in the above is that a divide by zero occurs for vertical part of the wellbore path. Hence we will use the trick to compute scale and convergence using a dummy survey (see tutorial) only for the first point and the last point, and then output as follows in a new property for v4 of the API

    "scaleConvergence": [
        {
            "scalefactor": 1.000241,
            "convergence": 0.51355,
            "point": {
                "x": 1999999.99999999,
                "y": 9999999.99999969,
                "z": 0.0
            }
        },
        {
            "scalefactor": 1.000243,
            "convergence": 0.51354,
            "point": {
                "x": 2000206.0812087534,
                "y": 9999956.440082304,
                "z": -1181.1945438868763
            }
        }
     ]

Ad 2)

This method is documented in section 3.3 of the wellbore doc which is linked in the [CRS Convert tutorial] (https://community.opengroup.org/osdu/platform/system/reference/crs-conversion-service/-/blob/master/docs/v3/tutorial/CRS_Convert_Service_howto.md#5-computing-a-wellbore-trajectory-from-directional-survey-data).

The "GNL" method requires the input survey observables to be grid north referenced, i.e.,

method: "GNL"
inputKind: "MD_Incl_Azim"
azimuthReference: "GridNorth" (or "GN" - however this is coded, probably detected by the first letter is g or G). (TBD if this is required. It may be possible that TN is given if the AzimuthalEquidistant method may work with it and a projected CRS).

It is a very simple method that notionally works as follows. First, as always, the minimum curvature computed local offsets are computed (which are "true to scale"). Then these are simply added to the 3D surface location in the same projected CRS. By not scaling these local offsets back to the map projection the difference between a cubical coordinates and curved coordinates are ignored.

However, it is desirable to output the scale factor, and that would not be possible if the minimum curvature offsets are simply added to the starting location. Hence, the implementation of GNL is to internally call "AzimuthalEquidistant", to compute scale factor, and "unscale" the results:

Check inputKind and azimuthReference and interpolate options.
If method: "GNL" then actually do "azimuthalEquidistant" internally.
Then calculate the psf at each station as per the above (ad 1).
Then "unscale" the XY trajectory as described here.

Acceptance criteria:

documentation in tutorial and api specification.
implementation (passed tests)

Ad 3)

Add optional input MD_i and output Stations_i
Implement case to deal with interpolation_interval = Number (e.g., 10)
Implemented, tested, accepted.
Example in tutorial

Minimum curvature interpolation is done at given MD_i. Math is described in section 2.3 of "OSDU_wellbore_calculations.docx" which is linked in the [CRS Convert tutorial] (https://community.opengroup.org/osdu/platform/system/reference/crs-conversion-service/-/blob/master/docs/v3/tutorial/CRS_Convert_Service_howto.md#5-computing-a-wellbore-trajectory-from-directional-survey-data).

The algorithm is summarized as:

First compute the minimum curvature offsets as normal with the stations that have MD,INC,AZI observables.
In a second pass, for each desired MD_i[i],
- a. Find the station before and after MD_i[i].
- b. Interpolate the Dog Leg with the equations provided at the desired MD_i.
- c. Compute the interpolated INC_i and AZI_i.
- d. Compute the local offsets dx,dy,dz.
- e. Add those offsets to the previous (real) station.
Output calculated (incl. interpolated) stations in an array stations_i.

To trigger interpolation at MD, require (check that):

inputKind=="MD_Incl_Azim"
MD_i[]: either
- a constant interpolation interval
- an array with "md_i" values as input with depths at which to interpolate the path.

Option 1 to interpolate the trajectory every 1 [unitZ] for MD_i as interval:

 {
    "MD_i": {
        "md_interval": 1.0
    }
}

Option 2 for specified MD_i values:

A bad request error with message should be thrown if both an interval and individual MD_i are in the request.

For directional survey:

  "inputStations": [
    {
      "md": 0,
      "inclination": 10,
      "azimuth": 100
    },
    {
      "md": 1000,
      "inclination": 20,
      "azimuth": 110
    },
    ...

The interpolated stations/densified path is returned as Stations_i as follows. (note the algorithm should set the "original" flag to true for MDs exactly at an inputStation, as shown for md=0 in the below example. For interpolated stations, the INC_i and AZI_i are returned as computed at the MD_i as shown).

    "stations_i": [
        {
            "md": 0.0,
            "inclination": 10.0,
            "azimuthTN": 100.51354989131318,
            "azimuthGN": 100.0,
            "dxTN": 0.0,
            "dyTN": 0.0,
            "scalefactor": 1.000241,
            "convergence": 0.51355,
            "point": {
                "x": 1999999.99999999,
                "y": 9999999.99999969,
                "z": 0.0
            },
            "wgs84Longitude": -85.88980921169528,
            "wgs84Latitude": 27.553258329196616,
            "dls": 0.0,
            "original": true,
            "dz": 0.0
        },
        {
            "md": 200.0,
            "inclination": 12.000031,
            "azimuthTN": 102.51354989131318,
            "azimuthGN": 102.0023012,
             ...
            "original": false,
...

Note:

Consider if we have MD, INC but not AZI if this interpolation can be used also to compute paths for INC-ONLY stations. In that algorithm one may want to interpolate the AZI.
Consider if this can be used to project out if MD_i is past the last MD (keep the same INC and AZI essentially as last station) in that special case.
Consider what to do if MD_i is given and interpolate==True. We propose just keep the algorithm to do whatever it does now and not change its behavior. One could consider this to flag to output interleave interpolated stations in stations.

Edited May 19, 2023 by KIRAN ALLAMSETY