NTv2 bicubic and bilinear interpolators are not currently 'reversible'; Add a (single) iteration to find grid shift values in the CRS of the grid

[jhaasdyk-au]

It is possible to identify locations in the NTv2 grid gda94_gda2020_conformal_and_distortion.gsb (at areas where grid shifts change quickly)
where the forward and reverse transformations are not reversible
when using GeodePy.transform.ntv2_2d() (either method = 'bilinear', or method = 'bicubic').
i.e. GDA94 => GDA2020 => GDA94 does not return the original coordinate to the expected precision.
This is caused by assuming that the input lat and long are in the CRS of the grid.

This could become problematic if repeated transformations are performed,
or if a process expects the original coordinate to be returned at a very high precision (e.g. if attempting to identifying nodes as identical when comparing datasets from different sources, e.g. if (re)providing a client with GDA94 data before an after a migration of primary datasets to GDA2020)

e.g. below using 'bilinear'; GDA94 => GDA2020 => GDA94, at
-31.6, 141.5 (in 'EAST' subgrid) => => -31.5999999992, 141.5000000011 (approx 0.15mm different per iteration)
-20.0, 122.6 (in 'WEST' subgrid => => -20.0000000000, 122.6000000083 (approx 0.87mm different per iteration)

This issue can be solved by iterating (when transforming in the reverse direction),
so that the grid shifts which are applied are derived at the nominal location in the CRS of the grid.
e.g. determine the shift at the GDA94 location, and apply those SAME shifts, whether in the forward OR reverse direction.

def ntv2_2d(...)
   ...
    if forward_tf:
        tf_lat = lat + shifts[0] / 3600
        tf_lon = lon - shifts[1] / 3600
    else:
        tf_lat = lat - shifts[0] / 3600
        tf_lon = lon + shifts[1] / 3600

      ##### ADD THE FOLLOWING ####
      # iterate to obtain shifts at the nominal coords in the CRS of the grid
      # e.g. gda94 location for gda94_gda2020_conformal_and_distortion.gsb
        shifts = interpolate_ntv2(ntv2_grid, tf_lat, tf_lon, method=method)
        if shifts[0] is None:
            raise ValueError('Coordinate outside of grid extents')
        tf_lat = lat - shifts[0] / 3600
        tf_lon = lon + shifts[1] / 3600