search

OSGeo / gdal

GDAL is an open source MIT licensed translator library for raster and vector geospatial data formats.
https://gdal.org
Other
4.76k stars 2.5k forks source link

gdalinfo crash on calculating statistics for Esri Tile Package (.tpkx) file #10229

Closed jgrocha closed 2 months ago

jgrocha commented 2 months ago

What is the bug?

Opening this *.tpkx file crashes QGIS with core dumped. The last log is:

src/core/providers/gdal/qgsgdalprovider.cpp:3091 : (bandStatistics) [0ms] theStats = 11 supportedStats = 110111
src/core/providers/gdal/qgsgdalprovider.cpp:3100 : (bandStatistics) [0ms] Using GDAL statistics.
src/core/providers/gdal/qgsgdalprovider.cpp:3116 : (bandStatistics) [0ms] bApproxOK = 1
src/core/providers/gdal/qgsgdalprovider.cpp:3129 : (bandStatistics) [0ms] myerval = 2
src/core/providers/gdal/qgsgdalprovider.cpp:3134 : (bandStatistics) [0ms] Calculating statistics by GDAL
Aborted (core dumped)

Running just gdalinfo -stats 88f936c4477a40a2.tpkx crashes.

Band 1 Block=256x256 Type=Byte, ColorInterp=Red
GDAL: Use hashset band block cache
Aborted (core dumped)

The file was exported from https://tiles.arcgis.com

image

I can reexport with other options. I just exported one level (18) to convert it to GeoTiff.

Steps to reproduce the issue

1) Download the example file from https://braga.geomaster.pt/nextcloud/index.php/s/TGoZfDGXeQSWcQ2/download/88f936c4477a40a2.tpkx 2) Run gdalinfo -stats 88f936c4477a40a2.tpkx

Versions and provenance

GDAL 3.9.0dev-f3432cf5e0, released 2024/04/06 (debug build)

Additional context

Besides the crash, there are things not working as expected, but I'm quite new to this format. This is the first time I'm processing a tpkx file. Sorry for any incorrect interpretation.

1) The bounding box is not properly detected by QGIS. 2) The reported size if quite huge Size is 67108863, 67108863. 3) I'm not able to export to an usable GeoTIFF file, passing the correct bounding box:

gdal_translate --config CHECK_DISK_FREE_SPACE NO -projwin -881694.0145031642 4480211.252062512 -803593.763394229 4438225.711865369 -of Gtiff -co TILED=YES -co "compress=LZW" -co NUM_THREADS=ALL_CPUS 88f936c4477a40a2.tpkx ortos06.tiff
jratike80 commented 2 months ago

No crash for me on Windows (GDAL 3.9.0, released 2024/05/07 from OSGeo4W). The image is rather large (67108863, 67108863) compared to the file size. Is it somehow very sparse or what?

gdalinfo 88f936c4477a40a2.tpkx -stats
ERROR 1: 88f936c4477a40a2.tpkx, band 1: Failed to compute statistics, no valid pixels found in sampling.
ERROR 1: 88f936c4477a40a2.tpkx, band 2: Failed to compute statistics, no valid pixels found in sampling.
ERROR 1: 88f936c4477a40a2.tpkx, band 3: Failed to compute statistics, no valid pixels found in sampling.
ERROR 1: 88f936c4477a40a2.tpkx, band 4: Failed to compute statistics, no valid pixels found in sampling.
Driver: ESRIC/Esri Compact Cache
Files: 88f936c4477a40a2.tpkx
Size is 67108863, 67108863
Coordinate System is:
PROJCRS["WGS 84 / Pseudo-Mercator",
    BASEGEOGCRS["WGS 84",
...
jratike80 commented 2 months ago

Gdalinfo does not find the image extent from the .tpkx. It seems to print the extent of the whole EPSG:3857. For GDAL the layer is an 8bit RGBA one.

According to ArcGIS layer properties the image is a single band int32 without a color map, but the image is shown with colors and the pixel query tool shows normal 8bit RGBA values. I think that there is some bug in the ArcGIS Pro Layer properties tool.

ArcGIS Pro shows also this:

The size is:
columns 75901
rows 63659

Statistics are missing.
Extent:
Top     4478165.415044 m
Bottom 4440150.533917 m
Left -865442.723857 m
Right -820117.357570 m

It may be that even this extent is mostly empty space and there are image data in much smaller area.

jgrocha commented 2 months ago

@jratike80 Thank you for testing on your side.

The first thing I noticed is the failing bounding box detection, by GDAL. It is a small image in Algarve, south of Portugal. The size and extent reported by ArcGIS Pro is correct.

Seems like the crash in Linux goes into glibc, far beyond my knowledge, like a thread related problem.

$ gdb gdalinfo
(gdb) run -stats 88f936c4477a40a2.tpkx
(...)
Lower Right (20037507.751,-20037507.751) (179d59'59.98"E, 85d 3' 4.06"S)
Center      (  -0.2958325,   0.2958325) (  0d 0' 0.01"W,  0d 0' 0.01"N)
Band 1 Block=256x256 Type=Byte, ColorInterp=Red

Program received signal SIGABRT, Aborted.
__pthread_kill_implementation (no_tid=0, signo=6, threadid=140737168106368) at ./nptl/pthread_kill.c:44
44      ./nptl/pthread_kill.c: No such file or directory.
(gdb)
rouault commented 2 months ago

I'm working on a fix. With it, when displaying the .tpkx in QGIS with a OSM background, I see a shift of about 100 m in easting/northing between both layers. I can't spot an obvious bug in the driver (especially since that shift doesn't match a whole tile). Is the original dataset properly georeferenced ? Does that shift also occurs in ArcGIS ?

jgrocha commented 2 months ago

I'm working on a fix. With it, when displaying the .tpkx in QGIS with a OSM background, I see a shift of about 100 m in easting/northing between both layers. I can't spot an obvious bug in the driver (especially since that shift doesn't match a whole tile). Is the original dataset properly georeferenced ? Does that shift also occurs in ArcGIS ?

Nice work @rouault. Glad to see it opening in QGIS :-)

You are right! The original image is not properly georeferenced. That's why I'm going back to 2006 to fix the image position.

jgrocha commented 2 months ago

Tested and working. Very well done @rouault 

gdalinfo 88f936c4477a40a2.tpkx
Driver: ESRIC/Esri Compact Cache
Files: 88f936c4477a40a2.tpkx
Size is 75901, 63659
Coordinate System is:
PROJCRS["WGS 84 / Pseudo-Mercator",
    BASEGEOGCRS["WGS 84",
        ENSEMBLE["World Geodetic System 1984 ensemble",
            MEMBER["World Geodetic System 1984 (Transit)"],
            MEMBER["World Geodetic System 1984 (G730)"],
            MEMBER["World Geodetic System 1984 (G873)"],
            MEMBER["World Geodetic System 1984 (G1150)"],
            MEMBER["World Geodetic System 1984 (G1674)"],
            MEMBER["World Geodetic System 1984 (G1762)"],
            MEMBER["World Geodetic System 1984 (G2139)"],
            MEMBER["World Geodetic System 1984 (G2296)"],
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ENSEMBLEACCURACY[2.0]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["Popular Visualisation Pseudo-Mercator",
        METHOD["Popular Visualisation Pseudo Mercator",
            ID["EPSG",1024]],
        PARAMETER["Latitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["easting (X)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["northing (Y)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Web mapping and visualisation."],
        AREA["World between 85.06°S and 85.06°N."],
        BBOX[-85.06,-180,85.06,180]],
    ID["EPSG",3857]]
Data axis to CRS axis mapping: 1,2
Origin = (-865442.723856601864100,4478166.012208217754960)
Pixel Size = (0.597164283559817,-0.597164283559817)
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  ( -865442.724, 4478166.012) (  7d46'27.86"W, 37d16'46.94"N)
Lower Left  ( -865442.724, 4440151.131) (  7d46'27.86"W, 37d 0'26.97"N)
Upper Right ( -820117.358, 4478166.012) (  7d22' 2.06"W, 37d16'46.94"N)
Lower Right ( -820117.358, 4440151.131) (  7d22' 2.06"W, 37d 0'26.97"N)
Center      ( -842780.041, 4459158.572) (  7d34'14.96"W, 37d 8'37.39"N)
Band 1 Block=256x256 Type=Byte, ColorInterp=Red
  Mask Flags: PER_DATASET ALPHA 
Band 2 Block=256x256 Type=Byte, ColorInterp=Green
  Mask Flags: PER_DATASET ALPHA 
Band 3 Block=256x256 Type=Byte, ColorInterp=Blue
  Mask Flags: PER_DATASET ALPHA 
Band 4 Block=256x256 Type=Byte, ColorInterp=Alpha

ESRI tpkx in QGIS

image