Error importing TIFF file using raster in R: "Failure during raster IO"

Question

I've been trying to import a TIFF formatted raster image into R using the raster::raster() function with no luck, and I would like to know if the file is corrupted, or if there is something I am missing here (e.g. GDAL or raster setting).

library(raster)
elevation <- raster('raster_image.tif')

A RasterLayer object is created successfully, but when I try to plot it I get this error.

> plot(elevation)
Error in rgdal::getRasterData(con, offset = offs, region.dim = reg, band = object@data@band) : 
Failure during raster IO

A link to download the .tif file is here: https://www.dropbox.com/s/fpxdzpi6ikngr0b/raster_image.tif?dl=0

Also, if I try and access its contents via indexing, I can get output for some values:

> elevation[1:3]
[1] 24.20100 24.25353 24.30607

but other elements I get the same error as for plot():

> elevation[1:300]
Error in rgdal::getRasterData(con, offset = offs, region.dim = reg, band = object@data@band) : 
  Failure during raster IO

• R Version: 4.0.3
• rgdal Version: 1.5-18
• GDAL Version: 3.0.4, released 2020/01/28

For context, the TIFF file was downloaded from the publicly available ELVIS elevation data repository in Australia. The actual TIFF is for a box around (151.0642° east, -33.9399° south) and is the 2020 version.

Things I have tried so far

I tried to use GDAL directly via the command line:

gdal_translate raster_image.tif output.tif                              

But I get the following error related to LZW compression:

Input file size is 2000, 2000                                             
0...10.ERROR 1: LZWDecode:Not enough data at scanline 1024 (short 240 bytes)   
ERROR 1: TIFFReadEncodedTile() failed.
ERROR 1: raster_image.tif, band 1: IReadBlock failed at X offset 2, Y offset 1: TIFFReadEncodedTile() failed.

I did manage to render it in QGIS where it doesn't appear to be corrupted:

but when I use the export tool to save to a GeoTIFF the resulting file has a large region missing.

Additionally, there is a warning in the "Rendering" tab of Log Messages (have removed actual file path):

2021-03-09T10:48:34     WARNING    raster_image_39b1418a_1419_4ea7_94ca_2c1dff6d2067 :: file/path/to/raster_image.tif, band 1: IReadBlock failed at X offset 2, Y offset 1: TIFFReadEncodedTile() failed.

I've also tried the gdalwarp solution in this related question:

>gdalwarp raster_image.tif raster_image-fix.tif
Creating output file that is 2000P x 2000L.
Processing raster_image.tif [1/1] : 0Using internal nodata values (e.g. -9999) for image raster_image.tif.
Copying nodata values from source raster_image.tif to destination raster_image-fix.tif.
ERROR 1: LZWDecode:Not enough data at scanline 1024 (short 240 bytes)
ERROR 1: TIFFReadEncodedTile() failed.
ERROR 1: raster_image.tif, band 1: IReadBlock failed at X offset 2, Y offset 1: TIFFReadEncodedTile() failed. 

but the resulting file when loaded into R has completely missing values:

> elevation2 <- raster('raster_image-fix.tif')
Warning message:
In showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj = prefer_proj) :
  Discarded datum Unknown based on GRS80 ellipsoid in CRS definition
> summary(elevation2)
        raster_image.fix
Min.                  NA
1st Qu.               NA
Median                NA
3rd Qu.               NA
Max.                  NA
NA's             3994240
Warning message:
In .local(object, ...) :
  summary is an estimate based on a sample of 1e+05 cells (2.5% of all cells)

Answer 1

Here's some python code that may help you or any other users with the same problem in the future. Using only builtin libs (no pyQGIS or arcpy) and GDAL_Translate binary - make sure that the GDAL install binary file location is in your %PATH%:

import os, sys, subprocess

baseFolder = 'A:\\'                                               # Change to match your base folder
tempFile   = os.path.join(os.environ.get('TEMP'),'delete_me.tif') # Temp convert file, overwritten on each iteration
fileCount  = 0

print('Process begins in {0:<80}'.format(baseFolder)),

# Open a null object to pipe the standard output and standard error
with open(os.devnull,'w') as OutNul:
    with open(os.devnull,'w') as ErrNul:
        
        # Open a text file for writing the file names of convert fails
        with open(r'A:\BadFiles.txt','w') as BadFileReport:
            
            # walk down from the top folder, this lists in every subfolder
            for (fullPath, Dirs, Files) in os.walk(baseFolder):
                
                # iterate individual files
                for thisFile in Files:
                    fileCount+= 1
                    
                    # decide if the file has a raster extension
                    if os.path.splitext(thisFile.lower())[1] in ['.tif','.img','.hdr']:
                        print '\rFile {0} {1:<80}'.format(fileCount,thisFile),
                        
                        # get the full path to this raster, fullPath from the
                        # walk() is the current folder being traversed
                        inRaster = os.path.join(fullPath,thisFile)
                        
                        # start GDAL_Translate in a subprocess, overwriting the temp raster and piping
                        # standard output and standard error to nothing
                        PObj     = subprocess.Popen(['GDAL_TRANSLATE','-of','GTIFF',inRaster,tempFile,
                                                                                 '-co', 'NUM_THREADS=ALL_CPUS', '-co', 'ZLEVEL=1',
                                                                                 '-CO', 'BIGTIFF=YES', '-co', 'compress=LZW'],
                                                                                 '-outsize','10%','10%', stdout=OutNul, stderr=ErrNul)
                        
                        PObj.wait() # wait until the process finishes
                        
                        if PObj.returncode != 0: # retuned object from the process has a non 0 if the process fails
                            # write this files' full path to the text file
                            BadFileReport.write('{}\n'.format(inRaster))
                            print('\r{0:<80} is bad'.format(thisFile))

# clean up the temporary raster
if os.path.exists(tempFile):
    os.remove(tempFile)

It's not quite crash-and-burn as calling subprocess.Popen and having that operation fail does not raise an exception, if the GDAL_Translate operation fails then that can be detected in the Popen object using returnCode, which is 0 for success and not 0 for failed operation. This operation will take a long time if you've got very big or lots of rasters, I've set the output to 10% of the original (10 x current cell size) to speed up the process and reduce free space requirement.