Fix identification of TAB file in GDA2020 CRS

Fixes #33007 (also requires GDAL master)
This commit is contained in:
Nyall Dawson 2019-11-28 11:54:01 +10:00
parent a4715a30c9
commit 3f94d55384
4 changed files with 80 additions and 38 deletions

View File

@ -784,8 +784,15 @@ bool QgsCoordinateReferenceSystem::createFromWkt( const QString &wkt )
#if PROJ_VERSION_MAJOR>=6
if ( d->mPj )
{
const QString authName( proj_get_id_auth_name( d->mPj.get(), 0 ) );
const QString authCode( proj_get_id_code( d->mPj.get(), 0 ) );
QString authName( proj_get_id_auth_name( d->mPj.get(), 0 ) );
QString authCode( proj_get_id_code( d->mPj.get(), 0 ) );
if ( authName.isEmpty() || authCode.isEmpty() )
{
// try 2, use proj's identify method and see if there's a nice (singular) candidate we can use
QgsProjUtils::identifyCrs( d->mPj.get(), authName, authCode );
}
if ( !authName.isEmpty() && !authCode.isEmpty() )
{
if ( loadFromAuthCode( authName, authCode ) )
@ -939,44 +946,18 @@ bool QgsCoordinateReferenceSystem::createFromProj4( const QString &proj4String )
if ( crs )
{
//crs = QgsProjUtils::crsToSingleCrs( crs.get() ) ;
int *confidence = nullptr;
if ( PJ_OBJ_LIST *crsList = proj_identify( QgsProjContext::get(), crs.get(), nullptr, nullptr, &confidence ) )
QString authName;
QString authCode;
if ( QgsProjUtils::identifyCrs( crs.get(), authName, authCode ) )
{
const int count = proj_list_get_count( crsList );
int bestConfidence = 0;
QgsProjUtils::proj_pj_unique_ptr matchedCrs;
for ( int i = 0; i < count; ++i )
const QString authid = QStringLiteral( "%1:%2" ).arg( authName, authCode );
if ( createFromOgcWmsCrs( authid ) )
{
if ( confidence[i] >= bestConfidence )
{
// prefer EPSG codes for compatibility with earlier qgis conversions
QgsProjUtils::proj_pj_unique_ptr candidateCrs( proj_list_get( QgsProjContext::get(), crsList, i ) );
candidateCrs = QgsProjUtils::crsToSingleCrs( candidateCrs.get() );
const QString authName( proj_get_id_auth_name( candidateCrs.get(), 0 ) );
if ( confidence[i] > bestConfidence || authName == QLatin1String( "EPSG" ) )
{
bestConfidence = confidence[i];
matchedCrs = std::move( candidateCrs );
}
}
}
proj_list_destroy( crsList );
proj_int_list_destroy( confidence );
if ( matchedCrs && bestConfidence >= 70 )
{
const QString authName( proj_get_id_auth_name( matchedCrs.get(), 0 ) );
const QString authCode( proj_get_id_code( matchedCrs.get(), 0 ) );
if ( !authName.isEmpty() && !authCode.isEmpty() )
{
const QString authid = QStringLiteral( "%1:%2" ).arg( authName, authCode );
if ( createFromOgcWmsCrs( authid ) )
{
locker.changeMode( QgsReadWriteLocker::Write );
if ( !sDisableProj4Cache )
sProj4Cache()->insert( proj4String, *this );
return true;
}
}
locker.changeMode( QgsReadWriteLocker::Write );
if ( !sDisableProj4Cache )
sProj4Cache()->insert( proj4String, *this );
return true;
}
}
}

View File

@ -186,6 +186,46 @@ QgsProjUtils::proj_pj_unique_ptr QgsProjUtils::crsToSingleCrs( const PJ *crs )
#endif
}
bool QgsProjUtils::identifyCrs( const PJ *crs, QString &authName, QString &authCode )
{
authName.clear();
authCode.clear();
if ( !crs )
return false;
int *confidence = nullptr;
if ( PJ_OBJ_LIST *crsList = proj_identify( QgsProjContext::get(), crs, nullptr, nullptr, &confidence ) )
{
const int count = proj_list_get_count( crsList );
int bestConfidence = 0;
QgsProjUtils::proj_pj_unique_ptr matchedCrs;
for ( int i = 0; i < count; ++i )
{
if ( confidence[i] >= bestConfidence )
{
// prefer EPSG codes for compatibility with earlier qgis conversions
QgsProjUtils::proj_pj_unique_ptr candidateCrs( proj_list_get( QgsProjContext::get(), crsList, i ) );
candidateCrs = QgsProjUtils::crsToSingleCrs( candidateCrs.get() );
const QString authName( proj_get_id_auth_name( candidateCrs.get(), 0 ) );
if ( confidence[i] > bestConfidence || authName == QLatin1String( "EPSG" ) )
{
bestConfidence = confidence[i];
matchedCrs = std::move( candidateCrs );
}
}
}
proj_list_destroy( crsList );
proj_int_list_destroy( confidence );
if ( matchedCrs && bestConfidence >= 70 )
{
authName = QString( proj_get_id_auth_name( matchedCrs.get(), 0 ) );
authCode = QString( proj_get_id_code( matchedCrs.get(), 0 ) );
}
}
return !authName.isEmpty() && !authCode.isEmpty();
}
bool QgsProjUtils::coordinateOperationIsAvailable( const QString &projDef )
{
if ( projDef.isEmpty() )

View File

@ -103,6 +103,14 @@ class CORE_EXPORT QgsProjUtils
*/
static proj_pj_unique_ptr crsToSingleCrs( const PJ *crs );
/**
* Attempts to identify a \a crs, matching it to a known authority and code within
* an acceptable level of tolerance.
*
* Returns TRUE if a matching authority and code was found.
*/
static bool identifyCrs( const PJ *crs, QString &authName, QString &authCode );
/**
* Returns TRUE if a coordinate operation (specified via proj string) is available.
*/

View File

@ -50,6 +50,7 @@ class TestQgsCoordinateReferenceSystem: public QObject
void createFromSrid();
void sridCache();
void createFromWkt();
void createFromWktWithIdentify();
void fromWkt();
void wktCache();
void createFromESRIWkt();
@ -342,6 +343,18 @@ void TestQgsCoordinateReferenceSystem::createFromWkt()
QCOMPARE( myCrs.srsid(), GEOCRS_ID );
}
void TestQgsCoordinateReferenceSystem::createFromWktWithIdentify()
{
#if PROJ_VERSION_MAJOR>=6
QgsCoordinateReferenceSystem crs;
// see https://github.com/qgis/QGIS/issues/33007, WKT string from a MapInfo tabfile in GDA2020 projection
crs.createFromWkt( QStringLiteral( "BOUNDCRS[SOURCECRS[PROJCRS[\"unnamed\",BASEGEOGCRS[\"unnamed\",DATUM[\"Geocentric Datum of Australia 2020\",ELLIPSOID[\"GRS 80\",6378137,298.257222101,LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]]],PRIMEM[\"Greenwich\",0,ANGLEUNIT[\"degree\",0.0174532925199433,ID[\"EPSG\",9122]]]],CONVERSION[\"UTM zone 55S\",METHOD[\"Transverse Mercator\",ID[\"EPSG\",9807]],PARAMETER[\"Latitude of natural origin\",0,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8801]],PARAMETER[\"Longitude of natural origin\",147,ANGLEUNIT[\"degree\",0.0174532925199433],ID[\"EPSG\",8802]],PARAMETER[\"Scale factor at natural origin\",0.9996,SCALEUNIT[\"unity\",1],ID[\"EPSG\",8805]],PARAMETER[\"False easting\",500000,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8806]],PARAMETER[\"False northing\",10000000,LENGTHUNIT[\"metre\",1],ID[\"EPSG\",8807]],ID[\"EPSG\",17055]],CS[Cartesian,2],AXIS[\"easting\",east,ORDER[1],LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]],AXIS[\"northing\",north,ORDER[2],LENGTHUNIT[\"metre\",1,ID[\"EPSG\",9001]]]]],TARGETCRS[GEOGCRS[\"WGS 84\",DATUM[\"World Geodetic System 1984\",ELLIPSOID[\"WGS 84\",6378137,298.257223563,LENGTHUNIT[\"metre\",1]]],PRIMEM[\"Greenwich\",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[\"unknown\"],AREA[\"World\"],BBOX[-90,-180,90,180]],ID[\"EPSG\",4326]]],ABRIDGEDTRANSFORMATION[\"Transformation to WGS84\",METHOD[\"Position Vector transformation (geog2D domain)\",ID[\"EPSG\",9606]],PARAMETER[\"X-axis translation\",-0.06155,ID[\"EPSG\",8605]],PARAMETER[\"Y-axis translation\",0.01087,ID[\"EPSG\",8606]],PARAMETER[\"Z-axis translation\",0.04019,ID[\"EPSG\",8607]],PARAMETER[\"X-axis rotation\",-0.0394924,ID[\"EPSG\",8608]],PARAMETER[\"Y-axis rotation\",-0.0327221,ID[\"EPSG\",8609]],PARAMETER[\"Z-axis rotation\",-0.0328979,ID[\"EPSG\",8610]],PARAMETER[\"Scale difference\",1.000000009994,ID[\"EPSG\",8611]]]]" ) );
QVERIFY( crs.isValid() );
QCOMPARE( crs.authid(), QStringLiteral( "EPSG:7855" ) );
#endif
}
void TestQgsCoordinateReferenceSystem::fromWkt()
{
QgsCoordinateReferenceSystem myCrs = QgsCoordinateReferenceSystem::fromWkt( geoWkt() );