QGIS show gdal_polygonized file in EPSG:54004 instead of EPSG:3857

Question

I am polygonizing a TIFF with the gdal_polygonize.py and the resulting Shape is shown in QGIS in the wrong projection - EPSG:54004 instead of EPSG:3857.

I have a TIFF in EPSG:3857 and polygonize it with

gdal_polygonize.py input.tif -f "ESRI Shapefile" output.shp

After doing this, the corresponding .prj file looks ok to me and I also checked it with various other sources like this

PROJCS["WGS_84_Pseudo_Mercator",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Mercator"],PARAMETER["central_meridian",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["Meter",1],PARAMETER["standard_parallel_1",0.0]]

When I load the output.shp into QGIS it automatically detects the file in EPSG:54004 which results in an offset compared with the source. After changing the Shapes projection in QGIS manually to EPSG:3857 both files, the TIFF and the Shape, do align.

Could this be a problem with QGIS or a problem with GDAL?

Offset between TIFF and Shape in the red rectangle

Update

gdalsrsinfo on the original TIFF shows

PROJ.4 : '+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs'

OGC WKT :
PROJCS["WGS 84 / Pseudo-Mercator",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs"],
    AUTHORITY["EPSG","3857"]]

Answer 1

It is the obscurity of Pseudo Mercator that leads to the offset. Both projections (3857 and 54004) share the same WKT definition, but it is treated differently.

Google (Pseudo) Mercator takes lat/lon coordinates of the ellipsoid, and uses them as they were on a sphere. Hence the different definitions a=b= 6378137 in the proj string vs SPHEROID["WGS 84",6378137,298.257223563 in WKT.

EPSG:3857 is the only projection method that handles it correct, but if no EPSG code is given, QGIS sometimes misinterprets the .prj file and assigns EPSG:54004, which leads to the offset you noticed.

There is no other way than assigning EPSG:3857 with Set CRS for Layer, then using Save As ... which adds a .qpj file, with the EPSG code number incorporated. This is a workaround for QGIS, hence it is not incorporated or read in GDAL tools like gdal_polygonize.py