search

qgis / QGIS

QGIS is a free, open source, cross platform (lin/win/mac) geographical information system (GIS)
https://qgis.org
GNU General Public License v2.0
10.59k stars 3k forks source link

QGIS 3.10.0-2 (proj 6.2.1) reprojects layers incorrectly while 3.10.0-1 (proj 5.2.0) does it okey. #33121

Closed denchat closed 4 years ago

denchat commented 4 years ago

Describe the bug

When add google XYZ tiles into a UTM project, the tiles appears to be over-shifted to the right to the UTM layers.

How to Reproduce

Use QGIS-OSGeo4W-3.10.0-2-Setup-x86_64.exe or download nightly OSGeo4W qgis-rel-dev. It seems both executables change their dependencies from PROJ 5.2.0 to PROJ 6.2.1.

QGIS and OS versions

QGIS:

OS:

gioman commented 4 years ago

@denchat duplicate of https://github.com/qgis/QGIS/issues/30569 ?

denchat commented 4 years ago

I'm not sure if this issue and #33129 are from the same cause. So I attached some files to be investigated further.

For my case, it appears that incorrect reprojections occur with on-the-fly transformation between WGS 1984 and Indian 1975 datum

This attachment _[magnolias_waterfront_epsg4240.zip]_ contains:

The project contain a landmark in Indian 1975 and google satellite tiles in WGS84 datum. Project CRS is in Indian 1975 datum.

When the project is opened in 3.10.0-1 [6c816b4204], the landmark is in its place. :

correct_reprojected

But when the project is opened in 3.10.0-2 [6ffa89eb3e] or nightly release, the landmark appears shifted. Like the WGS 1984 datum layer are on-the-fly reprojected incorrectly :

incorrect_reprojected

nyalldawson commented 4 years ago

Likely already fixed with 3.10.1

denchat commented 4 years ago

@nyalldawson I'm sorry for asking this. I could not find 3.10.1 in OSGeo4W yet. Even qgis-dev is 3.11.0-26.

Is there a way to obtain windows binary of 3.10.1?

gioman commented 4 years ago

Is there a way to obtain windows binary of 3.10.1?

is not released yet.

denchat commented 4 years ago

@gioman Could you just reopen this, please? It turns out that 3.10.1 is still affected.

nirvn commented 4 years ago

I can confirm the issue over here. In the console, those following critical messages are being printed again and again:

2019-12-08T12:29:52     CRITICAL    Cannot normalize transform between EPSG:3857 and EPSG:4240
2019-12-08T12:29:53     CRITICAL    Cannot normalize transform between EPSG:4240 and EPSG:3857
nirvn commented 4 years ago

One quick follow up that'd be good to do: test ogr2ogr to reproject a pseudomercador point coordinate into EPSG 4240 coordinate to rule whether the issue lies with QGIS or GDAL.

nyalldawson commented 4 years ago

@nirvn I don't see that -- are you self building proj? Maybe you need to update your build

nirvn commented 4 years ago

@nyalldawson , tested yesterday using osgeo4win's latest qgis-dev package.

dr1542 commented 4 years ago

I confirm this issue, Qgis versions begins from 3.10.0-2 (proj 6.2.1) transforms layers with too big difference. See project with test layers in attacment. Last running as expected QGIS version was 3.10.0-1 (proj 5.2.0)

crs_test.zip

nirvn commented 4 years ago

@dr1542 , have you confirmed that the regression is still happening in latest 3.10.1 released a few days ago?

dr1542 commented 4 years ago

Yes, of course, i tried today in

QGIS version 3.10.1-A Coruña QGIS code revision ef24c526da
Compiled against Qt 5.11.2 Running against Qt 5.11.2
Compiled against GDAL/OGR 3.0.2 Running against GDAL/OGR 3.0.2
Compiled against GEOS 3.8.0-CAPI-1.13.1 Running against GEOS 3.8.0-CAPI-1.13.1
Compiled against SQLite 3.29.0 Running against SQLite 3.29.0
PostgreSQL Client Version 11.5 SpatiaLite Version 4.3.0
QWT Version 6.1.3 QScintilla2 Version 2.10.8
Compiled against PROJ 6.2.1 Running against PROJ Rel. 6.2.1, November 1st, 2019
OS Version Windows 7 SP 1 (6.1)
dr1542 commented 4 years ago

And night build 3.11.0-Master (proj 6.3.0) runs wiht the same (wrong) result.

roya0045 commented 4 years ago

Does it only occur between WGS84 and Indian 1975?

denchat commented 4 years ago

Pyramide_du_Louvre_epsg23031.zip

@roya0045 I tried WGS84 and ED50 there. It seems not affected. It looks okey.

nirvn commented 4 years ago

Did some more looking into this issue today. If there is a bug, it's GDAL-related, not QGIS.

Long story short, GDAL 2 and GDAL 3 converts the following WGS84 point (100.510492, 13.727615) into two different values: image

@rouault , was GDAL 2 wrong all this time and GDAL 3 right, or is it the other way around? :)

dr1542 commented 4 years ago

This looks like a bug in proj: it doesn't seem to use the +to_wgs parameter in the CRS description. I tried to create two similar user CRS with the same parameters, but in one of them used to_wgs parameter, and applied them to two copies of vector layer. When i turned on reprojection to project crs (wgs84-based) the objects of these two layers was placed at the same position, but when i opened this project in older version (qgis3.10.0.1, proj5.0.2) objects is displaced according to to_wgs parameters. Maybe the experts here will find the right explanation for this case, but in my opinion it seems like a proj bug.

nirvn commented 4 years ago

On ubuntu, with proj 5.2, the following command returns the correct coordinates:

echo 100.510492 13.727615 0 0 | cct +proj=longlat +a=6377276.345 +b=6356075.41314024 +towgs84=210,814,289,0,0,0,0 +no_defs
100.5137742335   13.7259483134       41.5030        0.0000

Rebooting on Windows to check against proj 6.

nirvn commented 4 years ago

Under proj 6, using the old proj string (which isn't what is returned when triggering EPSG:4240), the result matches proj 5:

C:\>cct +proj=longlat +a=6377276.345 +b=6356075.41314024 +towgs84=210,814,289,0,0,0,0 +no_defs
100.510492 13.727615 0 0
100.5137742335   13.7259483134       41.5030        0.0000

However, is I use the new string for EPSG:4240 under proj 6, the result shifts:

C:\>cct +proj=longlat +ellps=evrst30 +towgs84=293,836,318,0.5,1.6,-2.8,2.1 +no_defs
100.510492 13.727615 0 0
100.5154435957   13.7254888019       14.8954        0.0000
dr1542 commented 4 years ago

Hmm, I tried this test: echo 37.55 55.6 0 0 | .\cct +proj=tmerc +lat_0=0 +lon_0=38.48333333333333 +k=1 +x_0=2250000 +y_0=-5712900.566 +ellps=krass +towgs84=24.47,-130.89,-81.56,0,0,0.13,-0.22 +units=m +no_defs for proj5.2.0 and proj6.3.0

with the identical correct result: 2191283.0094 451622.7530 -7.0342 0.0000

So I concluded that there is no bug in proj.

nirvn commented 4 years ago

@dr1542 , let's try not to spam the issue here. The original author raised a reprojection issue with EPSG:4240.

dr1542 commented 4 years ago

Sorry, I think that there is more common case than just reprojection between EPSG:4240 and EPSG:3857. Other similar issues: #30569, #32795, #32928

nirvn commented 4 years ago

@dr1542 , by "similar" you refer to distinct projection-related issues. While those all fall within a same projection theme, it's better to keep those distinct as the underlying issues/bugs could be entirely different.

denchat commented 4 years ago

@nirvn There 4 conversion parameter sets of "toWGS84" on Indian 1975. In affected QGIS version, even though you change one conversion parameter to another, it seems there is nothing happens.

Perhaps, it could be about how the affected version treats the conversion setting in project file properties. The thing is why some conversions are affected while some (e.g. ED50) are not!?

Another thing, it looks like PROJ v6 migrates into SQLite for all its data store. Nevertheless, the helmert_transformation_table, which contain "toWGS84" parameter sets, looks all intact. So I guess it could be more likely there might be something with the client side than with the PROJ itself.

dr1542 commented 4 years ago

I checked coordinate transformations inside qgis for gdal python api and qgis api for old and new versions qgis between WGS84 and user crs with +towgs84 parameter and without it by running simple test script in qgis python console. Results for old qgis, proj, gdal versions : Gdal (version 2040100) transforms WGS84 point(37.55 55.6) to USER crs with +towgs84 POINT (2191283.00944974 451622.752971655 -7.03415696974844) to USER crs without +towgs84 POINT (2191166.68962418 451630.580835653 0) QGIS 3.10.0-A Coruña (6c816b4204) transforms WGS84 point(37.55 55.6) to USER crs with +towgs84 Point (2191283.0094497287645936 451622.75297106057405472) to USER crs without +towgs84 Point (2191166.68962418008595705 451630.58083565160632133) There is no difference between GDAL and QGIS transformations Results for new versions: Gdal (version 3010000) transforms WGS84 point(37.55 55.6) to USER crs with +towgs84 POINT (2191283.00944973 451622.752968665 0) to USER crs without +towgs84 POINT (2191166.68962418 451630.580833254 0) QGIS 3.11.0-Master (2b195e9450) transforms WGS84 point(37.55 55.6) to USER crs with +towgs84 Point (2191166.68962418427690864 451630.58083325438201427) to USER crs without +towgs84 Point (2191166.68962418427690864 451630.58083325438201427)

A completely different result! Two last lines shows than new qgis not used towgs84 parameter in transformation and no suspicion for GDAL.

Then i tried next simple test for EPSG4240 proj4 definition string as other sample with towgs84 parameter. In python console I typed: print(QgsCoordinateReferenceSystem().fromProj4('+proj=longlat +ellps=evrst30 +towgs84=293,836,318,0.5,1.6,-2.8,2.1 +no_defs').toWkt()) QGIS 3.10.0-A Coruña (6c816b4204) returns GEOGCS["Everest 1830",DATUM["unknown",SPHEROID["evrst30",6377276.345,300.8017],TOWGS84[293,836,318,0.5,1.6,-2.8,2.1]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]] The same definition in WKT But QGIS 3.11.0-Master (2b195e9450): GEOGCS["Indian 1975",DATUM["Indian_1975",SPHEROID["Everest 1830 (1937 Adjustment)",6377276.345,300.8017,AUTHORITY["EPSG","7015"]],AUTHORITY["EPSG","6240"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4240"]] Was returned WKT definition with missed TOWGS84 parameter!

So it looks like qgis losses towgs84 parameter when converts crs to internal form.

Test script for gdal and qgis transformations in attachment. transformTest.zip

nirvn commented 4 years ago

@dr1542 , great work coming up with this script.

nirvn commented 4 years ago

Some more "fun", along the lines of what @dr1542 did in his useful post.

Using gdal's python API, running the following:

crs = osr.SpatialReference()
crs.ImportFromEPSG(4240)
crs.ExportToWkt()

will output:

'GEOGCS["Indian 1975",DATUM["Indian_1975",SPHEROID["Everest 1830 (1937 Adjustment)",6377276.345,300.8017,AUTHORITY["EPSG","7015"]],TOWGS84[293,836,318,0.5,1.6,-2.8,2.1],AUTHORITY["EPSG","6240"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4240"]]'

Running the following QGIS python API:

crs = QgsCoordinateReferenceSystem('EPSG:4240')
crs.toWkt()

returns:

'GEOGCS["Indian 1975",DATUM["Indian_1975",SPHEROID["Everest 1830 (1937 Adjustment)",6377276.345,300.8017,AUTHORITY["EPSG","7015"]],AUTHORITY["EPSG","6240"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4240"]]'

As far as I can see, both are by default outputting PJ_WKT1_GDAL strings, yet as @dr1542 pointed out, the towgs84 parameter is gone missing in QGIS' returned wkt string.

nyalldawson commented 4 years ago

@nirvn

Check projinfo "EPSG:4240":

WKT2:2019 string:
GEOGCRS["Indian 1975",
    DATUM["Indian 1975",
        ELLIPSOID["Everest 1830 (1937 Adjustment)",6377276.345,300.8017,
            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["Thailand - onshore and Gulf of Thailand"],
        BBOX[5.63,97.34,20.46,105.64]],
    ID["EPSG",4240]]

I can't explain gdal's behavior, but we get the result directly from proj...

nyalldawson commented 4 years ago

Although hmm, there's a difference there...

nyalldawson commented 4 years ago

@nirvn how did you get your results? (what build, etc)

I see:

crs = QgsCoordinateReferenceSystem('EPSG:4240')
crs.toWkt(QgsCoordinateReferenceSystem.WKT2_2018)

gives

'GEOGCRS["Indian 1975",DATUM["Indian 1975",ELLIPSOID["Everest 1830 (1937 Adjustment)",6377276.345,300.8017]],CS[ellipsoidal,2],AXIS["geodetic latitude (Lat)",north],AXIS["geodetic longitude (Lon)",east],UNIT["degree",0.0174532925199433],USAGE[SCOPE["unknown"],AREA["Thailand - onshore and Gulf of Thailand"],BBOX[5.63,97.34,20.46,105.64]],ID["EPSG",4240]]'

i.e. identical to projinfo

nirvn commented 4 years ago

@nyalldawson , I'm on windows, using osgeo4w's latest qgis-dev package.

On my machine, if I run projinfo -o WKT1_GDAL "EPSG:4240" , I get:

WKT1_GDAL:
GEOGCS["Indian 1975",
    DATUM["Indian_1975",
        SPHEROID["Everest 1830 (1937 Adjustment)",6377276.345,300.8017,
            AUTHORITY["EPSG","7015"]],
        TOWGS84[293,836,318,0.5,1.6,-2.8,2.1],
        AUTHORITY["EPSG","6240"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4240"]]

The parameters match that of QGIS, except for TOWGS84 being absent.

nyalldawson commented 4 years ago

https://github.com/OSGeo/PROJ/issues/1794 relates

denchat commented 4 years ago

Not sure that this is related with the issue

QgsCoordinateTransformContext https://qgis.org/pyqgis/master/core/QgsCoordinateTransformContext.html

addSourceDestinationDatumTransform(self, sourceCrs: QgsCoordinateReferenceSystem, destinationCrs: QgsCoordinateReferenceSystem, sourceTransformId: int, destinationTransformId: int) → bool

Deprecated since version Has: no effect on builds based on Proj 6.0 or later, use addCoordinateOperation() instead.

Is it possible that some how the call graph invokes addSourceDestinationDatumTransform() instead of addCoordinateOperation() which it should since 3.10.0-2 after adopting PROJ 6?

denchat commented 4 years ago

@nyalldawson I has another strange case here

I try avoiding the issue of inconsistant tranform of 4240. So I create a new custom CRS for each of them using their own proj string.

I reproject the landmark 13.727615N 100.510492E to those custom CRSes. That is

I assign project CRS to EPSG:4326 WGS84, then save and exit QGIS (3.10.0-1)

Now I try opening the project with different QGIS versions.

  1. When open with 3.10.0-1 PROJ 5.2
    • custom crs transform 1304: QGIS interprets it as EPSG:4240 (as expected because 1812 doesn't exist yet)
    • custom crs transform 1812: QGIS interprets it as a custom crs.

All 3 points share the exact location.

proj5

  1. When opened with 3.10.1-1 PROJ 6.2
    • custom crs transform 1304: QGIS interprets it as a custom crs. (as expected because 1812 now already exists)
    • custom crs transform 1812: QGIS interprets it as a EPSG:4240

However instead of sharing the exact location, the transform 1304 point is shifted to the right and down!?

proj6

But when I change project CRS to ESPG:4240, the shifted tranform 1304 point moves back to the left and up like we saw before.

proj6-4240

_magnolias_waterfrontepsg4240.zip

I'm confused that how come changing from WGS84 to Indian 1975 can shift "tranform 1304" to the opposite side?

Please investigate this case also.

nyalldawson commented 4 years ago

@denchat should all be resolved with current master version -- please retest

denchat commented 4 years ago

@nyalldawson Confirm fixed on 3.11.0-40. This is such a marvellous work. You save me from the shift!

I'm now not seeing the backport needed tag. Is this already merged in 3.10-dev also?

nyalldawson commented 4 years ago

Yes, it was backported to 3.10 this morning

nyalldawson commented 4 years ago

And thanks for the confirmation and testing, it's very much appreciated!

nirvn commented 4 years ago

You save me from the shift!

Heh.