2

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: tiff rendering in QGIS

but when I use the export tool to save to a GeoTIFF the resulting file has a large region missing. exported image not properly rendering

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)
4
  • 5
    It looks like the raster is corrupt, I downloaded it and copied with ArcGIS to an ERDAS IMG format and there's a huge chunk missing. The block size (GDALInfo) is 512 x 512, one block is bad. Note the CRS is MGA2020 Zone 56 not geographic. The reason it renders in QGIS is because the GeoTIFF file has internal overviews (GeoTIFF can have external or internal pyramids), when zoomed to full extent you're only seeing an overview, when you zoom in the data will be missing. I suggest you download the data again and if that doesn't fix it contact the custodian of the data. Mar 9 at 5:38
  • Thank you @MichaelStimson for your response. This explains everything that I've seen, and yes I got some weird effects zooming in to the image in QGIS. I tried downloading the data again but I had the exact same issue, so I've contacted Geoscience Australia to let them know. Out of interest, do you know of an open-source tool to identify missing data blocks in TIF files, similar to how you did it with ArcGIS? Mar 10 at 2:09
  • Do you have any python ability? Calling GDAL_Translate in subprocess.Popen then polling the process for its return code would identify the bad ones - GDAL_Translate return 0 for a good file and 1 for a bad file. If you're interested I could answer with some working code. Mar 10 at 3:08
  • Thanks, I ended up doing a similar thing using the try() function in R, where try(expr) returns an object of class "try-error" if the the expr gives an error. I am interested to see the python code, so if you don't mind writing it up as an answer that would be great. Mar 10 at 23:04
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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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