Skip to content

Commit 90a9f3f

Browse files
rouaultrouault
authored
Merge pull request OSGeo#3729 from rouault/crs_to_crs_geocentric_crs
cs2cs / proj_create_crs_to_crs(): fix regression with geocentric CRS
2 parents 852f3c4 + 190bfb2 commit 90a9f3f

File tree

5 files changed

+306
-55
lines changed

5 files changed

+306
-55
lines changed

‎src/4D_api.cpp

Lines changed: 194 additions & 51 deletions
Original file line numberDiff line numberDiff line change
@@ -245,8 +245,24 @@ int pj_get_suggested_operation(PJ_CONTEXT *,
245245
const auto &alt = opList[i];
246246
bool spatialCriterionOK = false;
247247
if (direction == PJ_FWD) {
248-
if (coord.xyzt.x >= alt.minxSrc && coord.xyzt.y >= alt.minySrc &&
249-
coord.xyzt.x <= alt.maxxSrc && coord.xyzt.y <= alt.maxySrc) {
248+
if (alt.pjSrcGeocentricToLonLat) {
249+
if (alt.minxSrc == -180 && alt.minySrc == -90 &&
250+
alt.maxxSrc == 180 && alt.maxySrc == 90) {
251+
spatialCriterionOK = true;
252+
} else {
253+
PJ_COORD tmp = coord;
254+
pj_fwd4d(tmp, alt.pjSrcGeocentricToLonLat);
255+
if (tmp.xyzt.x >= alt.minxSrc &&
256+
tmp.xyzt.y >= alt.minySrc &&
257+
tmp.xyzt.x <= alt.maxxSrc &&
258+
tmp.xyzt.y <= alt.maxySrc) {
259+
spatialCriterionOK = true;
260+
}
261+
}
262+
} else if (coord.xyzt.x >= alt.minxSrc &&
263+
coord.xyzt.y >= alt.minySrc &&
264+
coord.xyzt.x <= alt.maxxSrc &&
265+
coord.xyzt.y <= alt.maxySrc) {
250266
spatialCriterionOK = true;
251267
} else if (alt.srcIsLonLatDegree && coord.xyzt.y >= alt.minySrc &&
252268
coord.xyzt.y <= alt.maxySrc) {
@@ -264,8 +280,24 @@ int pj_get_suggested_operation(PJ_CONTEXT *,
264280
}
265281
}
266282
} else {
267-
if (coord.xyzt.x >= alt.minxDst && coord.xyzt.y >= alt.minyDst &&
268-
coord.xyzt.x <= alt.maxxDst && coord.xyzt.y <= alt.maxyDst) {
283+
if (alt.pjDstGeocentricToLonLat) {
284+
if (alt.minxDst == -180 && alt.minyDst == -90 &&
285+
alt.maxxDst == 180 && alt.maxyDst == 90) {
286+
spatialCriterionOK = true;
287+
} else {
288+
PJ_COORD tmp = coord;
289+
pj_fwd4d(tmp, alt.pjDstGeocentricToLonLat);
290+
if (tmp.xyzt.x >= alt.minxDst &&
291+
tmp.xyzt.y >= alt.minyDst &&
292+
tmp.xyzt.x <= alt.maxxDst &&
293+
tmp.xyzt.y <= alt.maxyDst) {
294+
spatialCriterionOK = true;
295+
}
296+
}
297+
} else if (coord.xyzt.x >= alt.minxDst &&
298+
coord.xyzt.y >= alt.minyDst &&
299+
coord.xyzt.x <= alt.maxxDst &&
300+
coord.xyzt.y <= alt.maxyDst) {
269301
spatialCriterionOK = true;
270302
} else if (alt.dstIsLonLatDegree && coord.xyzt.y >= alt.minyDst &&
271303
coord.xyzt.y <= alt.maxyDst) {
@@ -320,6 +352,14 @@ int pj_get_suggested_operation(PJ_CONTEXT *,
320352
}
321353

322354
//! @cond Doxygen_Suppress
355+
/**************************************************************************************/
356+
PJCoordOperation::~PJCoordOperation() {
357+
/**************************************************************************************/
358+
proj_destroy(pj);
359+
proj_destroy(pjSrcGeocentricToLonLat);
360+
proj_destroy(pjDstGeocentricToLonLat);
361+
}
362+
323363
/**************************************************************************************/
324364
bool PJCoordOperation::isInstantiable() const {
325365
/**************************************************************************************/
@@ -1659,11 +1699,11 @@ static void reproject_bbox(PJ *pjGeogToCrs, double west_lon, double south_lat,
16591699
}
16601700

16611701
/*****************************************************************************/
1662-
static PJ *add_coord_op_to_list(int idxInOriginalList, PJ *op, double west_lon,
1663-
double south_lat, double east_lon,
1664-
double north_lat, PJ *pjGeogToSrc,
1665-
PJ *pjGeogToDst, bool isOffshore,
1666-
std::vector<PJCoordOperation> &altCoordOps) {
1702+
static PJ *add_coord_op_to_list(
1703+
int idxInOriginalList, PJ *op, double west_lon, double south_lat,
1704+
double east_lon, double north_lat, PJ *pjGeogToSrc, PJ *pjGeogToDst,
1705+
const PJ *pjSrcGeocentricToLonLat, const PJ *pjDstGeocentricToLonLat,
1706+
bool isOffshore, std::vector<PJCoordOperation> &altCoordOps) {
16671707
/*****************************************************************************/
16681708

16691709
double minxSrc;
@@ -1675,19 +1715,35 @@ static PJ *add_coord_op_to_list(int idxInOriginalList, PJ *op, double west_lon,
16751715
double maxxDst;
16761716
double maxyDst;
16771717

1678-
reproject_bbox(pjGeogToSrc, west_lon, south_lat, east_lon, north_lat,
1679-
minxSrc, minySrc, maxxSrc, maxySrc);
1680-
reproject_bbox(pjGeogToDst, west_lon, south_lat, east_lon, north_lat,
1681-
minxDst, minyDst, maxxDst, maxyDst);
1718+
if (pjSrcGeocentricToLonLat) {
1719+
minxSrc = west_lon;
1720+
minySrc = south_lat;
1721+
maxxSrc = east_lon;
1722+
maxySrc = north_lat;
1723+
} else {
1724+
reproject_bbox(pjGeogToSrc, west_lon, south_lat, east_lon, north_lat,
1725+
minxSrc, minySrc, maxxSrc, maxySrc);
1726+
}
1727+
1728+
if (pjDstGeocentricToLonLat) {
1729+
minxDst = west_lon;
1730+
minyDst = south_lat;
1731+
maxxDst = east_lon;
1732+
maxyDst = north_lat;
1733+
} else {
1734+
reproject_bbox(pjGeogToDst, west_lon, south_lat, east_lon, north_lat,
1735+
minxDst, minyDst, maxxDst, maxyDst);
1736+
}
16821737

16831738
if (minxSrc <= maxxSrc && minxDst <= maxxDst) {
16841739
const char *c_name = proj_get_name(op);
16851740
std::string name(c_name ? c_name : "");
16861741

16871742
const double accuracy = proj_coordoperation_get_accuracy(op->ctx, op);
1688-
altCoordOps.emplace_back(idxInOriginalList, minxSrc, minySrc, maxxSrc,
1689-
maxySrc, minxDst, minyDst, maxxDst, maxyDst,
1690-
op, name, accuracy, isOffshore);
1743+
altCoordOps.emplace_back(
1744+
idxInOriginalList, minxSrc, minySrc, maxxSrc, maxySrc, minxDst,
1745+
minyDst, maxxDst, maxyDst, op, name, accuracy, isOffshore,
1746+
pjSrcGeocentricToLonLat, pjDstGeocentricToLonLat);
16911747
op = nullptr;
16921748
}
16931749
return op;
@@ -1798,6 +1854,57 @@ static PJ *create_operation_to_geog_crs(PJ_CONTEXT *ctx, const PJ *crs) {
17981854
return opGeogToCrs;
17991855
}
18001856

1857+
/*****************************************************************************/
1858+
static PJ *create_operation_geocentric_crs_to_geog_crs(PJ_CONTEXT *ctx,
1859+
const PJ *geocentric_crs)
1860+
/*****************************************************************************/
1861+
{
1862+
assert(proj_get_type(geocentric_crs) == PJ_TYPE_GEOCENTRIC_CRS);
1863+
1864+
auto datum = proj_crs_get_datum_forced(ctx, geocentric_crs);
1865+
assert(datum);
1866+
auto cs = proj_create_ellipsoidal_2D_cs(ctx, PJ_ELLPS2D_LONGITUDE_LATITUDE,
1867+
nullptr, 0);
1868+
auto ellps = proj_get_ellipsoid(ctx, datum);
1869+
proj_destroy(datum);
1870+
double semi_major_metre = 0;
1871+
double inv_flattening = 0;
1872+
proj_ellipsoid_get_parameters(ctx, ellps, &semi_major_metre, nullptr,
1873+
nullptr, &inv_flattening);
1874+
// It is critical to set the prime meridian to 0
1875+
auto lon_lat_crs = proj_create_geographic_crs(
1876+
ctx, "unnamed crs", "unnamed datum", proj_get_name(ellps),
1877+
semi_major_metre, inv_flattening, "Reference prime meridian", 0,
1878+
nullptr, 0, cs);
1879+
proj_destroy(ellps);
1880+
proj_destroy(cs);
1881+
1882+
// Create the transformation from this geocentric CRS to the lon-lat one
1883+
auto operation_ctx = proj_create_operation_factory_context(ctx, nullptr);
1884+
proj_operation_factory_context_set_spatial_criterion(
1885+
ctx, operation_ctx, PROJ_SPATIAL_CRITERION_PARTIAL_INTERSECTION);
1886+
proj_operation_factory_context_set_grid_availability_use(
1887+
ctx, operation_ctx,
1888+
PROJ_GRID_AVAILABILITY_DISCARD_OPERATION_IF_MISSING_GRID);
1889+
auto op_list =
1890+
proj_create_operations(ctx, geocentric_crs, lon_lat_crs, operation_ctx);
1891+
proj_operation_factory_context_destroy(operation_ctx);
1892+
proj_destroy(lon_lat_crs);
1893+
1894+
const int nOpCount = op_list == nullptr ? 0 : proj_list_get_count(op_list);
1895+
if (nOpCount != 1) {
1896+
proj_context_log_debug(ctx, "Cannot compute transformation from "
1897+
"geocentric CRS to geographic CRS");
1898+
proj_list_destroy(op_list);
1899+
return nullptr;
1900+
}
1901+
1902+
auto op = proj_list_get(ctx, op_list, 0);
1903+
assert(op);
1904+
proj_list_destroy(op_list);
1905+
return op;
1906+
}
1907+
18011908
/*****************************************************************************/
18021909
PJ *proj_create_crs_to_crs(PJ_CONTEXT *ctx, const char *source_crs,
18031910
const char *target_crs, PJ_AREA *area) {
@@ -1846,21 +1953,46 @@ pj_create_prepared_operations(PJ_CONTEXT *ctx, const PJ *source_crs,
18461953
const PJ *target_crs, PJ_OBJ_LIST *op_list)
18471954
/*****************************************************************************/
18481955
{
1849-
auto pjGeogToSrc = create_operation_to_geog_crs(ctx, source_crs);
1850-
if (!pjGeogToSrc) {
1851-
proj_context_log_debug(ctx,
1852-
"Cannot create transformation from geographic "
1853-
"CRS of source CRS to source CRS");
1854-
return {};
1956+
PJ *pjGeogToSrc = nullptr;
1957+
PJ *pjSrcGeocentricToLonLat = nullptr;
1958+
if (proj_get_type(source_crs) == PJ_TYPE_GEOCENTRIC_CRS) {
1959+
pjSrcGeocentricToLonLat =
1960+
create_operation_geocentric_crs_to_geog_crs(ctx, source_crs);
1961+
if (!pjSrcGeocentricToLonLat) {
1962+
// shouldn't happen
1963+
return {};
1964+
}
1965+
} else {
1966+
pjGeogToSrc = create_operation_to_geog_crs(ctx, source_crs);
1967+
if (!pjGeogToSrc) {
1968+
proj_context_log_debug(
1969+
ctx, "Cannot create transformation from geographic "
1970+
"CRS of source CRS to source CRS");
1971+
return {};
1972+
}
18551973
}
18561974

1857-
auto pjGeogToDst = create_operation_to_geog_crs(ctx, target_crs);
1858-
if (!pjGeogToDst) {
1859-
proj_context_log_debug(ctx,
1860-
"Cannot create transformation from geographic "
1861-
"CRS of target CRS to target CRS");
1862-
proj_destroy(pjGeogToSrc);
1863-
return {};
1975+
PJ *pjGeogToDst = nullptr;
1976+
PJ *pjDstGeocentricToLonLat = nullptr;
1977+
if (proj_get_type(target_crs) == PJ_TYPE_GEOCENTRIC_CRS) {
1978+
pjDstGeocentricToLonLat =
1979+
create_operation_geocentric_crs_to_geog_crs(ctx, target_crs);
1980+
if (!pjDstGeocentricToLonLat) {
1981+
// shouldn't happen
1982+
proj_destroy(pjSrcGeocentricToLonLat);
1983+
proj_destroy(pjGeogToSrc);
1984+
return {};
1985+
}
1986+
} else {
1987+
pjGeogToDst = create_operation_to_geog_crs(ctx, target_crs);
1988+
if (!pjGeogToDst) {
1989+
proj_context_log_debug(
1990+
ctx, "Cannot create transformation from geographic "
1991+
"CRS of target CRS to target CRS");
1992+
proj_destroy(pjSrcGeocentricToLonLat);
1993+
proj_destroy(pjGeogToSrc);
1994+
return {};
1995+
}
18641996
}
18651997

18661998
try {
@@ -1888,18 +2020,21 @@ pj_create_prepared_operations(PJ_CONTEXT *ctx, const PJ *source_crs,
18882020

18892021
const bool isOffshore = areaName && strstr(areaName, "- offshore");
18902022
if (west_lon <= east_lon) {
1891-
op = add_coord_op_to_list(i, op, west_lon, south_lat, east_lon,
1892-
north_lat, pjGeogToSrc, pjGeogToDst,
1893-
isOffshore, preparedOpList);
2023+
op = add_coord_op_to_list(
2024+
i, op, west_lon, south_lat, east_lon, north_lat,
2025+
pjGeogToSrc, pjGeogToDst, pjSrcGeocentricToLonLat,
2026+
pjDstGeocentricToLonLat, isOffshore, preparedOpList);
18942027
} else {
18952028
auto op_clone = proj_clone(ctx, op);
18962029

1897-
op = add_coord_op_to_list(i, op, west_lon, south_lat, 180,
1898-
north_lat, pjGeogToSrc, pjGeogToDst,
1899-
isOffshore, preparedOpList);
2030+
op = add_coord_op_to_list(
2031+
i, op, west_lon, south_lat, 180, north_lat, pjGeogToSrc,
2032+
pjGeogToDst, pjSrcGeocentricToLonLat,
2033+
pjDstGeocentricToLonLat, isOffshore, preparedOpList);
19002034
op_clone = add_coord_op_to_list(
19012035
i, op_clone, -180, south_lat, east_lon, north_lat,
1902-
pjGeogToSrc, pjGeogToDst, isOffshore, preparedOpList);
2036+
pjGeogToSrc, pjGeogToDst, pjSrcGeocentricToLonLat,
2037+
pjDstGeocentricToLonLat, isOffshore, preparedOpList);
19032038
proj_destroy(op_clone);
19042039
}
19052040

@@ -1908,10 +2043,14 @@ pj_create_prepared_operations(PJ_CONTEXT *ctx, const PJ *source_crs,
19082043

19092044
proj_destroy(pjGeogToSrc);
19102045
proj_destroy(pjGeogToDst);
2046+
proj_destroy(pjSrcGeocentricToLonLat);
2047+
proj_destroy(pjDstGeocentricToLonLat);
19112048
return preparedOpList;
19122049
} catch (const std::exception &) {
19132050
proj_destroy(pjGeogToSrc);
19142051
proj_destroy(pjGeogToDst);
2052+
proj_destroy(pjSrcGeocentricToLonLat);
2053+
proj_destroy(pjDstGeocentricToLonLat);
19152054
return {};
19162055
}
19172056
}
@@ -2070,13 +2209,10 @@ PJ *proj_create_crs_to_crs_from_pj(PJ_CONTEXT *ctx, const PJ *source_crs,
20702209
}
20712210

20722211
const auto backup_errno = proj_context_errno(ctx);
2073-
if ((P == nullptr ||
2074-
(op_count == 1 &&
2075-
(!mayNeedToReRunWithDiscardMissing ||
2076-
errorIfBestTransformationNotAvailable ||
2077-
singleOpIsInstanciable == static_cast<int>(true))) ||
2078-
proj_get_type(source_crs) == PJ_TYPE_GEOCENTRIC_CRS ||
2079-
proj_get_type(target_crs) == PJ_TYPE_GEOCENTRIC_CRS)) {
2212+
if (P == nullptr ||
2213+
(op_count == 1 && (!mayNeedToReRunWithDiscardMissing ||
2214+
errorIfBestTransformationNotAvailable ||
2215+
singleOpIsInstanciable == static_cast<int>(true)))) {
20802216
proj_list_destroy(op_list);
20812217
ctx->forceOver = false;
20822218

@@ -2835,21 +2971,28 @@ static bool isSpecialCaseForWGS84_to_GDA2020(const std::string &opName) {
28352971
return opName.find("GDA2020 to WGS 84 (2)") != std::string::npos;
28362972
}
28372973

2838-
PJCoordOperation::PJCoordOperation(int idxInOriginalListIn, double minxSrcIn,
2839-
double minySrcIn, double maxxSrcIn,
2840-
double maxySrcIn, double minxDstIn,
2841-
double minyDstIn, double maxxDstIn,
2842-
double maxyDstIn, PJ *pjIn,
2843-
const std::string &nameIn, double accuracyIn,
2844-
bool isOffshoreIn)
2974+
PJCoordOperation::PJCoordOperation(
2975+
int idxInOriginalListIn, double minxSrcIn, double minySrcIn,
2976+
double maxxSrcIn, double maxySrcIn, double minxDstIn, double minyDstIn,
2977+
double maxxDstIn, double maxyDstIn, PJ *pjIn, const std::string &nameIn,
2978+
double accuracyIn, bool isOffshoreIn, const PJ *pjSrcGeocentricToLonLatIn,
2979+
const PJ *pjDstGeocentricToLonLatIn)
28452980
: idxInOriginalList(idxInOriginalListIn), minxSrc(minxSrcIn),
28462981
minySrc(minySrcIn), maxxSrc(maxxSrcIn), maxySrc(maxySrcIn),
28472982
minxDst(minxDstIn), minyDst(minyDstIn), maxxDst(maxxDstIn),
28482983
maxyDst(maxyDstIn), pj(pjIn), name(nameIn), accuracy(accuracyIn),
28492984
isOffshore(isOffshoreIn),
28502985
isPriorityOp(isSpecialCaseForNAD83_to_NAD83HARN(name) ||
28512986
isSpecialCaseForGDA94_to_WGS84(name) ||
2852-
isSpecialCaseForWGS84_to_GDA2020(name)) {
2987+
isSpecialCaseForWGS84_to_GDA2020(name)),
2988+
pjSrcGeocentricToLonLat(pjSrcGeocentricToLonLatIn
2989+
? proj_clone(pjSrcGeocentricToLonLatIn->ctx,
2990+
pjSrcGeocentricToLonLatIn)
2991+
: nullptr),
2992+
pjDstGeocentricToLonLat(pjDstGeocentricToLonLatIn
2993+
? proj_clone(pjDstGeocentricToLonLatIn->ctx,
2994+
pjDstGeocentricToLonLatIn)
2995+
: nullptr) {
28532996
const auto IsLonLatOrLatLon = [](const PJ *crs, bool &isLonLatDegreeOut,
28542997
bool &isLatLonDegreeOut) {
28552998
const auto eType = proj_get_type(crs);

‎src/iso19111/c_api.cpp

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -9110,7 +9110,8 @@ PJ *proj_normalize_for_visualization(PJ_CONTEXT *ctx, const PJ *obj) {
91109110
alt.idxInOriginalList, minxSrc, minySrc, maxxSrc,
91119111
maxySrc, minxDst, minyDst, maxxDst, maxyDst,
91129112
pjNormalized, co->nameStr(), alt.accuracy,
9113-
alt.isOffshore);
9113+
alt.isOffshore, alt.pjSrcGeocentricToLonLat,
9114+
alt.pjDstGeocentricToLonLat);
91149115
}
91159116
}
91169117
return pjNew.release();

0 commit comments

Comments
 (0)