0

When comparing in QGIS the output of gdal2tiles.py to the raster from which the tiles were generated there is a clear difference.

In this image Pink is the gdal2tiles.py output, green is the original raster: enter image description here

The raster is projected to EPSG:4326, the gdal2tiles command is: gdal2tiles.py {file} {output_dir} --s_srs 'EPSG:4326' --zoom=8-18 --processes={cores} --exclude --xyz This image from QGIS also has the project CRS in EPSG:4326.

Question:

I just want to know why this could be and what I can do so that the two sources line up?


I have looked at tools like gdalwarp, but it doesnt seem to work: gdalwarp -s_srs 'EPSG:4326' -tr 29404 36007 -tap input.tiff output.tiff

gdalinfo on the raster file returns:

Driver: GTiff/GeoTIFF
Files: wgs84_sb_pluvial_30y_demo.tiff
Size is 29404, 36007
Coordinate System is:
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]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (-2.030342225189838,51.716732655211317)
Pixel Size = (0.000051188440515,-0.000051188440515)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_SOFTWARE=GRASS GIS 7.8.5 with GDAL 3.1.4
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  -2.0303422,  51.7167327) (  2d 1'49.23"W, 51d43' 0.24"N)
Lower Left  (  -2.0303422,  49.8735905) (  2d 1'49.23"W, 49d52'24.93"N)
Upper Right (  -0.5251973,  51.7167327) (  0d31'30.71"W, 51d43' 0.24"N)
Lower Right (  -0.5251973,  49.8735905) (  0d31'30.71"W, 49d52'24.93"N)
Center      (  -1.2777698,  50.7951616) (  1d16'39.97"W, 50d47'42.58"N)
2
  • What do you try to reach with your gdalwarp command? -tr 29404 36007 would give pixels which are 29404 units wide and 36007 units high. I guess that -ts is what you mean to use gdal.org/programs/gdalwarp.html.
    – user30184
    Aug 8, 2022 at 17:10
  • Ah yes that was a mistake. I thought this may help while reprojecting between rasters between CRS EPSG:27700 to EPSG:4326, reducing the stretched effect
    – sobmortin
    Aug 8, 2022 at 17:57

1 Answer 1

1

If the issue that you mean is that the pixels are not aligned, that happens because gdal2tiles by default creates tiles to follow the Google compatible tiling schema https://gdal.org/programs/gdal2tiles.html

**-p <PROFILE>, --profile=<PROFILE>**

Tile cutting profile (mercator, geodetic, raster) - default ‘mercator’ (Google Maps compatible).

In the web Mercator tiling schema the limits of each tile and also each pixel are strictly defined beforehand. See images in https://docs.microsoft.com/en-us/bingmaps/articles/bing-maps-tile-system. The pixels of the source image and tiled target image can only match accurately if the extents of the source image by chance match exactly with the tiles. If they do not match then tiling will lead into resampling and the pixels will have a slight shift.

There should be no shift between the pixels at the zoom level of the base tiles if you use the tiling profile "raster" but then you won't have the same interoperability than with web Mercator tiles.

Excuse me if I understood wrong what is the issue.

5
  • No, I don't think you have misunderstood. I tried this: gdal2tiles.py input.tif output_dir --s_srs 'EPSG:4326' --zoom=8-18 --processes=3 --exclude --xyz --profile=raster but got the error: IndexError: list index out of range. I do need to display the output of gdal2tiles.py in mapbox, if that effects the choice of profile.
    – sobmortin
    Aug 8, 2022 at 17:53
  • It seems this error is because I can't run gdal2tiles.py with a greater zoom level than 8 when I run it with --profle=raster. Any idea why that would be?
    – sobmortin
    Aug 9, 2022 at 7:07
  • I guess that at zoom level 8 in the local raster profile the pixel size is the same as in the original image and it does not make sense to continue further. If original pixel size is 1 m splitting the pixels into 0.5 m, 0.25 m etc. pixels would only waste disk space.
    – user30184
    Aug 9, 2022 at 7:23
  • First, I think your original answer is correct in identifying the cause of the inconsistencies. Using the --profile=raster arg prevents me zooming further than 8 (I need to go to 18 in order to maintain resolution in mapbox. Picking up on your comment: The pixels of the source image and tiled target image can only match accurately if the extents of the source image by chance match exactly, is it possibel to ensure the source image exactly match with the tiled target?
    – sobmortin
    Aug 10, 2022 at 14:33
  • About being not able to get zoom levels which require oversampling with the raster profile, see github.com/OSGeo/gdal/issues/6207.
    – user30184
    Aug 17, 2022 at 6:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.