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

从 GeodePy.transform.ntv2_2d() 开始，跟踪其对 interpolate_ntv2() 的调用，并使用 gda94_gda2020_conformal_and_distortion.gsb 以及 EAST 和 WEST 示例作为参考案例。检查逆变换路径，并验证一次迭代是否能在名义网格 CRS 位置获取偏移量，然后确认 GDA94 → GDA2020 → GDA94 是否能以预期精度将坐标还原。

It is possible to identify locations in the NTv2 grid gda94_gda2020_conformal_and_distortion.gsb
where the forward and reverse transformations are not perfectly reversible in the current GeodePy NTv2 functions
i.e. 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 the current functions applying the input lat and long as if they are in the CRS of the grid.
However, where grid shift values change quickly between neighbouring nodes, even a 1-2m difference in inputs has a noticeable effect.

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 their 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