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

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

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