[GIS] How to write BIGTIFF raster files with rasterio

bigtiffpythonrasterio

A few months ago I had a hard time in finding a way to read BIGTIFF files. In the end I only managed to do so by creating a new conda environment with the following options:

conda create -n env_name -c conda-forge gdal vs2015_runtime=14

I must admit I have no idea what that instruction does (vs2015_runtime=14???), but it did the job. Now I am simply trying to merge a number of tif tiles into one unique tif raster file. I am using the code below:

...
...

#Merge tiles
mosaic, out_trans = merge(tiles)

#Update metadata
out_meta.update({"driver": "GTiff",
                 "height": mosaic.shape[1],
                 "width": mosaic.shape[2],
                 "nodata": -1,
                 "transform": out_trans,
                 "compress": 'lzw'
                })

#Output path
out_path = '...'

#Build mosaic
with rasterio.open(out_path, "w", **out_meta) as dest:
    dest.write(mosaic)
dest.close()

This code has actually worked fine on other occasions when trying to merge tif tiles, but this time it's giving me the following error messages:

ERROR 1: .\GB_OUTPUT\ALR\GB_COMBINED_ALR.tif, band 1: An error occurred while writing a dirty block from FlushCache
ERROR 1: TIFFAppendToStrip:Maximum TIFF file size exceeded. Use BIGTIFF=YES creation option.
ERROR 1: TIFFAppendToStrip:Maximum TIFF file size exceeded. Use BIGTIFF=YES creation option.
ERROR 1: TIFFAppendToStrip:Maximum TIFF file size exceeded. Use BIGTIFF=YES creation option.
ERROR 1: TIFFAppendToStrip:Maximum TIFF file size exceeded. Use BIGTIFF=YES creation option.

There are 2 errors there. One BIGTIFF related, and another one about some dirty block. Any idea why I can't create BIGTIFF files and how can I overcome the issue? Does the dirty block just mean there's a corrupt tif tile?

Best Answer

I tried with the following and it worked fine:

with rasterio.open(out_path, 
                   "w", 
                   **out_meta, 
                   tiled=True, 
                   blockxsize=256, 
                   blockysize=256, 
                   BIGTIFF='YES') as dest:

Related Solutions

[GIS] Reading raster files by block with rasterio

Reading Rasters by block can be done in rasterio and I'd argue it's easier than in GDAL. There is even a tutorial on windowed read/write over at GitHub.

Let's take a look at the read functions arguments, which allows you to set a window to read data from:

def read(self, indexes=None, out=None, window=None, masked=False,
        boundless=False):
"""Read raster bands as a multidimensional array
    Parameters
    ----------
    window : a pair (tuple) of pairs of ints, optional
        The optional `window` argument is a 2 item tuple. The first
        item is a tuple containing the indexes of the rows at which
        the window starts and stops and the second is a tuple
        containing the indexes of the columns at which the window
        starts and stops. For example, ((0, 2), (0, 2)) defines
        a 2x2 window at the upper left of the raster dataset.
"""

So all you need now is a way to create the window tuples. That's what rasterios block_windows function is for:

def block_windows(self, bidx=0):
    """Returns an iterator over a band's block windows and their
    indexes.
    [....]
    The primary use of this function is to obtain windows to pass to
    read_band() for highly efficient access to raster block data.
    """

So in your case I guess the code would look something like this:

with rasterio.open("your/data/geo.tif") as src:
     for block_index, window in src.block_windows(1):
         block_array = src.read(window=window)
         result_block = some_calculation(block_array)

[GIS] Open a raster file with read and write permissions in rasterio

In rasterio <= 1.0.7 this is not possible. You have to write and close, then read. If you want to write again, just overwrite the memfile object.

For example:

from rasterio.io import MemoryFile
import rasterio
import numpy as np

def create_memory_file(data, west_bound, north_bound, cellsize, driver='GTIFF'):
    #data is a numpy array
    if data.ndim ==2: # Handle 2 or 3D input arrays
        data = np.expand_dims(data, axis=0)
    dtype = data.dtype
    shape = data.shape
    transform = rasterio.transform.from_origin(west_bound, north_bound, cellsize, cellsize)
    with MemoryFile() as memfile:
        with memfile.open(
                driver=driver, width= shape[2], height = shape[1],
                transform=transform, count=shape[0], dtype=dtype) as dataset:
            dataset.write(data)
        return memfile.open()  # <==== Re-open the memfile


data = np.array([[1,2,3], [4,5,6], [7,8,9]]).astype(np.int32)
memfile = create_memory_file(data, 0, 2, 0.5)
print(memfile.read())

data = np.array([[4,5,6], [7,8,9], [10,11, 12]]).astype(np.int32)
memfile = create_memory_file(data, 0, 2, 0.5)
print(memfile.read())

Output:

[[[1 2 3]
  [4 5 6]
  [7 8 9]]]

[[[ 4  5  6]
  [ 7  8  9]
  [10 11 12]]]

As of rasterio 1.0.8 you can write and read as datasets in a MemoryFile are opened in r+/w+ mode. Note that the underlying GDAL driver needs to support creation of a new dataset from scratch not just creation of a copy of an existing dataset.