This also fixes conversion between geocentric latlong and geodetic latlong
with cs2cs. This was dealt with in PR 1093, but in the wrong direction
(the geocentric latitude must be <= in absolute value to the geodetic one)
The issue here was linked to the semantics of the +geoc specifier, which
affects the semantics of the input coordinates in the forward direction
(+geoc means that the input coordinate is is a geocentric latitude),
which defeats the logic of doing A to B by using the inverse path of A
and the forward path of B.

Relates to https://github.com/OSGeo/gdal/issues/3998

Before that change, cs2cs on a NAD83(HARN) to WGS84 transformation would
use the "NAD83(HARN) to WGS 84 (1)" transformation (a null Helmert
shift) that is valid for whole US, including non-CONUS areas, even when
used on points located on CONUS that has a "NAD83(HARN) to WGS 84 (3)"
transformation (non-null Helmert shift) with same accuracy (1m).

But if doing  EPSG:2874 "NAD83(HARN) / California zone 5 (ftUS)" to
WGS84, we would use this later "NAD83(HARN) to WGS 84 (3)"
transformation because the area of use of EPSG:2874 restricts to CONUS.
This isn't consistant.

With that change, we now have more consistent behavior, even if it can
be argued which of the 2 transformations is the best...

$ echo 34 -120 | src/cs2cs -d 8 EPSG:4326 "NAD83(HARN)" | src/cs2cs "NAD83(HARN)"  EPSG:2874
5955507.74	1828410.98 0.00

$ echo 34 -120 | src/cs2cs  EPSG:4326   EPSG:2874
5955507.74	1828410.98 0.00

If using the proposed patch in PR #2616 refreshed as,

```
diff --git a/src/4D_api.cpp b/src/4D_api.cpp
index c7725d3a..3bd58a97 100644
--- a/src/4D_api.cpp
+++ b/src/4D_api.cpp
@@ -998,7 +998,7 @@ static void reproject_bbox(PJ* pjGeogToCrs,
         constexpr int XY_SIZE = N_STEPS_P1 * 4;
         std::vector<double> x(XY_SIZE);
         std::vector<double> y(XY_SIZE);
-        const double step_lon = (east_lon - west_lon) / N_STEPS;
+        const double step_lon = (((east_lon >= west_lon) ? east_lon : east_lon + 360) - west_lon) / N_STEPS;
         const double step_lat = (north_lat - south_lat) / N_STEPS;
         for( int j = 0; j <= N_STEPS; j++ )
         {
@@ -1244,7 +1244,7 @@ std::vector<PJCoordOperation> pj_create_prepared_operations(PJ_CONTEXT *ctx,
             {
                 const bool isOffshore =
                     areaName && strstr(areaName, "- offshore");
-                if( west_lon <= east_lon )
+                if( true || west_lon <= east_lon )
                 {
                     op = add_coord_op_to_list(i, op,
                         west_lon, south_lat, east_lon, north_lat,
```

we'd get an incorrect result on one point, using the 'Pulkovo 1942 to WGS 84 (16)' less accurate transformation
due to mis-computation of bounding boxes.
```
--- ../test/cli/tv_out.dist	2021-04-07 22:29:02.706095239 +0200
+++ test/cli/tv_out	2021-04-07 22:29:35.809579495 +0200
@@ -499,7 +499,7 @@
 ##############################################################
 Check that we can use a transformation spanning the antimeridian (should use Pulkovo 1942 to WGS 84 (20))
 50 179.999999999	49d59'59.36"N	179d59'52.133"W 0.000
-50 -179.999999999	49d59'59.36"N	179d59'52.133"W 0.000
+50 -179.999999999	49d59'59.422"N	179d59'52.184"W 0.000
 ##############################################################
 Check that we can use a transformation spanning the antimeridian (should use Pulkovo 1942 to WGS 84 (20))
 5540944.47  499999.999	49d59'59.36"N	179d59'52.133"W 0.000
```