Is this exception “osgeo::proj::io::FactoryException” a bug? in 9.4.0

loftyscholar commented 1 month ago

When I compile and run the qgis project, I got this exception. Acturally, the proj string is "+proj=longlat +a=10000 +no_defs +type=crs"" When debuging source code, I found this exception is triggered in the function: std::string PROJStringParser::Private::guessBodyName(double a) {

auto ret = Ellipsoid::guessBodyName(dbContext_, a);
if (ret == "Non-Earth body" && dbContext_ == nullptr && ctx_ != nullptr) {
    dbContext_ =
        ctx_->get_cpp_context()->getDatabaseContext().as_nullable();
    if (dbContext_) {
        ret = Ellipsoid::guessBodyName(dbContext_, a);
    }
}
return ret;

} When calling, variable a is 10000.000000000000.

Actually, another function was called, i.e.,

std::string AuthorityFactory::identifyBodyFromSemiMajorAxis(double semi_major_axis, double tolerance) const { auto res = d->run("SELECT name, (ABS(semi_major_axis - ?) / semi_major_axis ) " "AS rel_error FROM celestial_body WHERE rel_error <= ?", {semi_major_axis, tolerance}); if (res.empty()) { throw FactoryException("no match found"); } if (res.size() > 1) { for (const auto &row : res) { if (row[0] != res.front()[0]) { throw FactoryException("more than one match found"); } } } return res.front()[0]; } When calling, res.size() is 3, and for the second row, FactoryException was thrown.

Example of problem

Demonstrate the problem you have found, either by using PROJ tools like cs2cs, cct or projinfo or add a code snippet that highlights the problem using the PROJ C/C++ API.

A "Minimal, Complete and Verifiable Example" will make it much easier for maintainers to help you: http://matthewrocklin.com/blog/work/2018/02/28/minimal-bug-reports

# Your code here

Problem description

[this should explain why the current behavior is a problem and why the expected output is a better solution.]

Expected Output

Environment Information

PROJ version (proj)
Operation System Information

Installation method

conda, apt-get, from source, etc...

rouault commented 1 month ago

When I compile and run the qgis project, I got this exception.

This exception might be emitted by PROJ with such an unsual semi-major axis, but it should be caught, at least in the C function that QGIS uses: which is it? With only the content of this report, it is not clear if there's actually something to change in PROJ or in QGIS Might be worth to open it rather as a QGIS issue with the description of the related QGIS scenario

loftyscholar commented 1 month ago

Thanks for your reply. This unsual semi-major axis is located in ellipsoid table in proj.db, and its code is "107911" and its name is "Ananke_2000_IAU_IAG".

rouault commented 1 month ago

@loftyscholar ok, I see indeed that's an issue when trying to do something like projinfo -s "+proj=longlat +a=10000 +no_defs +type=crs" -t ESRI:104911 which should be a no-op operation, because there are several celestial bodies that match the 10000 radius. But couldn't you just avoid using a deprecated PROJ.4 string like "+proj=longlat +a=10000 +no_defs +type=crs" and use the ESRI:104911 code in QGIS. That should avoid the issue.

rouault commented 1 month ago

or use the WKT2 representation:

GEOGCRS["GCS_Ananke_2000",
    DATUM["D_Ananke_2000",
        ELLIPSOID["Ananke_2000_IAU_IAG",10000,0,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Reference_Meridian",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Not known."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["ESRI",104911]]

Alternatively you may define the PROJ_IGNORE_CELESTIAL_BODY=YES environment variable to avoid an error like

bin/projinfo -s  "+proj=longlat +R=10000 +no_defs +type=crs" -t ESRI:104911
createOperations() failed with: Source and target ellipsoid do not belong to the same celestial body (Non-Earth body vs Ananke). You may override this check by setting the PROJ_IGNORE_CELESTIAL_BODY environment variable to YES.

loftyscholar commented 1 month ago

deprecated Ok, I'll have a try soon.

rouault commented 1 month ago

Improvements in PROJ queued in https://github.com/OSGeo/PROJ/pull/4149

With those, projinfo -s "+proj=longlat +R=10000 +no_defs +type=crs" -t ESRI:104911 now works out of the box

loftyscholar commented 1 month ago

Improvements in PROJ queued in #4149

With those, projinfo -s "+proj=longlat +R=10000 +no_defs +type=crs" -t ESRI:104911 now works out of the box

Thanks for your reply, have all the changes been commited to the proj source code? Where should I get the changes?

rouault commented 1 month ago

have all the changes been commited to the proj source code?

This has just been committed in the master branch, and will be soon in the 9.4 branch as well. This will be in the 9.4.1 release

OSGeo / PROJ