gdalwarp to change block size also changes pixel values

Question

I'm working with a GeoTIFF from http://www.fao.org/geonetwork/srv/en/main.home?uuid=ba4526fd-cdbf-4028-a1bd-5a559c4bff38

The original has data values from 0-11:

$ gdalinfo -hist ./data/FAO/glc_shv10_DOM.Tif
Driver: GTiff/GeoTIFF
Size is 43200, 21600
Metadata:
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_SOFTWARE=Arc/Info
  TIFFTAG_XRESOLUTION=119.999
  TIFFTAG_YRESOLUTION=119.999
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (    0.0,    0.0)
Lower Left  (    0.0,21600.0)
Upper Right (43200.0,    0.0)
Lower Right (43200.0,21600.0)
Center      (21600.0,10800.0)
Band 1 Block=43200x3 Type=Byte, ColorInterp=Gray
  256 buckets from -0.5 to 255.5:
  181470400 977689 29454097 33129685 71015243 18337497 3857498 122159 20544229
  28053726 11242330 534915447 0 0 0 ... (more zeros)

I tried to reorganize it with a 256x256 block size for subsequent steps in my processing pipeline, using:
gdalwarp -ot Byte -co COMPRESS=DEFLATE -co TILED=YES -s_srs "EPSG:4326" -t_srs "EPSG:4326" -to SRC_METHOD=NO_GEOTRANSFORM glc_shv10_DOM.Tif new_file.tif

I had used "-to SRC_METHOD=NO_GEOTRANSFORM" because the original file lacks information about its coordinate system but looks like EPSG:4326.

The new file has very different pixel values: 0, 19, 26, 30, 34, 37, 39, 41, 43, 45, 47, 48.

Driver: GTiff/GeoTIFF
Size is 43200, 21600
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]],
    USAGE[
        SCOPE["unknown"],
        AREA["World"],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
Origin = (0.000000000000000,21600.000000000000000)
Pixel Size = (1.000000000000000,-1.000000000000000)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_SOFTWARE=Arc/Info
  TIFFTAG_XRESOLUTION=119.999
  TIFFTAG_YRESOLUTION=119.999
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (       0.000,   21600.000) (  0d 0' 0.01"E,Invalid angle)
Lower Left  (   0.0000000,   0.0000000) (  0d 0' 0.01"E,  0d 0' 0.01"N)
Upper Right (   43200.000,   21600.000) (Invalid angle,Invalid angle)
Lower Right (   43200.000,       0.000) (Invalid angle,  0d 0' 0.01"N)
Center      (   21600.000,   10800.000) (Invalid angle,Invalid angle)
Band 1 Block=256x256 Type=Byte, ColorInterp=Gray
  256 buckets from -0.5 to 255.5:
  181470400 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 977689 0 0 0 0 0 0
  29454097 0 0 0 33129685 0 0 0 71015243 0 0 18337497 0 3857498 0
  122159 0 20544229 0 28053726 0 11242330 534915447 0 0 (more zeros)

I don't understand why the pixel values changed. Other questions where data values changed involve a scale_factor in the metadata, but my file does not have a scale_factor.

I also need to figure out what I'm missing about adding a coordinate system, the Corner Coordinates are still pixels but now being interpreted as degrees.

Answer 1

I downloaded the glc_shv10_DOM.Tif and get very different gdalinfo output with Origin = (-180.000000000000000,90.000000001440014) and Pixel Size = (0.008333333333400,-0.008333333333400). You appear to have stripped the georeferencing out of the tif somehow or possibly an issue with your GDAL install.

If you just want to tile internally and specify the CRS, I'd use gdal_translate with -a_srs to specify the existing CRS and -co TILED=YES to tile using GeoTIFF creation options (256x256 is the default, to specify something else use the block size options -co BLOCKXSIZE=n -co BLOCKYSIZE=n).

In the example below I gdal_translate directly from the downloaded zip file (using /vsizip/, one of GDALs virtual file systems). Note that the histogram and statistics of the input and output are the same.

# Tile it with gdal_translate, assume WGS84 (EPSG:4326), compress with LZW to match orig tif
$ gdal_translate -a_srs EPSG:4326 -co COMPRESS=LZW -co TILED=YES /vsizip/GlcShare_v10_Dominant.zip/glc_shv10_DOM.Tif glc_shv10_DOM.Tif

# Input
$ gdalinfo -hist -stats /vsizip/GlcShare_v10_Dominant.zip/glc_shv10_DOM.Tif

    Driver: GTiff/GeoTIFF
    Files: /vsizip/GlcShare_v10_Dominant.zip/glc_shv10_DOM.Tif
           /vsizip/GlcShare_v10_Dominant.zip/glc_shv10_DOM.aux
    Size is 43200, 21600
    Coordinate System is:
    LOCAL_CS["unnamed",
        UNIT["unknown",1]]
    Origin = (-180.000000000000000,90.000000001440014)
    Pixel Size = (0.008333333333400,-0.008333333333400)
    Metadata:
      </snip...>
    Image Structure Metadata:
      COMPRESSION=LZW
      INTERLEAVE=BAND
    Corner Coordinates:
    Upper Left  (-180.0000000,  90.0000000) 
    Lower Left  (-180.0000000, -90.0000000) 
    Upper Right ( 180.0000000,  90.0000000) 
    Lower Right ( 180.0000000, -90.0000000) 
    Center      (   0.0000000,   0.0000000) 
    Band 1 Block=43200x1 Type=Byte, ColorInterp=Gray
      Description = Band_1
      Min=0.000 Max=11.000 
      Minimum=0.000, Maximum=11.000, Mean=7.472, StdDev=4.571
      256 buckets from -0.5 to 255.5:
      181470400 977689 29454097 33129685 71015243 18337497 3857498 122159 20544229 28053726 11242330 534915447 </snip...>
      Metadata:
          </snip...>

# Output
$ gdalinfo -hist -stats glc_shv10_DOM.Tif

    Driver: GTiff/GeoTIFF
    Files: glc_shv10_DOM.Tif
    Size is 43200, 21600
    Coordinate System is:
    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"]]
    Origin = (-180.000000000000000,90.000000001440014)
    Pixel Size = (0.008333333333400,-0.008333333333400)
    Metadata:
      </snip...>
    Image Structure Metadata:
      COMPRESSION=LZW
      INTERLEAVE=BAND
    Corner Coordinates:
    Upper Left  (-180.0000000,  90.0000000) (180d 0' 0.00"W, 90d 0' 0.00"N)
    Lower Left  (-180.0000000, -90.0000000) (180d 0' 0.00"W, 90d 0' 0.00"S)
    Upper Right ( 180.0000000,  90.0000000) (180d 0' 0.00"E, 90d 0' 0.00"N)
    Lower Right ( 180.0000000, -90.0000000) (180d 0' 0.00"E, 90d 0' 0.00"S)
    Center      (   0.0000000,   0.0000000) (  0d 0' 0.00"E,  0d 0' 0.00"N)
    Band 1 Block=256x256 Type=Byte, ColorInterp=Gray
      Description = Band_1
      Min=0.000 Max=11.000 
      Minimum=0.000, Maximum=11.000, Mean=7.472, StdDev=4.571
      256 buckets from -0.5 to 255.5:
      181470400 977689 29454097 33129685 71015243 18337497 3857498 122159 20544229 28053726 11242330 534915447 </snip...>
      Metadata:
          </snip...>