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

pythonrasterio

I am trying to write a function which will help me construct a rasterio memory file and return the dataset reader object so I can then play around with it as needed. I am trying to avoid with statements because I would like to allow the user to read and write as they deem necessary. This is what I have so far:

from rasterio import mask, features, warp
from rasterio.io import MemoryFile
import rasterio
import numpy as np

def create_memory_file(data, west_bound, north_bound, cellsize, driver='AAIGrid'):
    #data is a numpy array
    dtype = data.dtype
    shape = data.shape
    transform = rasterio.transform.from_origin(west_bound, north_bound, cellsize, cellsize)
    memfile = MemoryFile()
    dataset = memfile.open(driver=driver, width= shape[1], height = shape[0], transform=transform, count=1, dtype = dtype)
    dataset.write(data,1)

    return dataset
def close_memory_file(memfile):
    memfile.close()


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

The last line throws an error:

Traceback (most recent call last):
  File "D:/11202750-002_RA2CE/Basis/common.py", line 214, in <module>
    memfile.read(1)
  File "rasterio\_io.pyx", line 209, in rasterio._io.DatasetReaderBase.read
rasterio.errors.UnsupportedOperation: not readable

The memfile apparently only has write permissions because of the following:

>>> memfile
<open BufferedDatasetWriter name='/vsimem/3c0572dd-a295-4f2b-a65d-22a404dd5d0f.' mode='w'>
>>> memfile.mode
'w'

I have tried various combinations of adding mode=r+, permissions=r+ to my memfile.open() call but it won't allow this.

How can I add the read/write mode at the level of the open statement?

Best Answer

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.

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:
        dataset = memfile.open(
                driver=driver, width= shape[2], height = shape[1],
                transform=transform, count=shape[0], dtype=dtype)
        dataset.write(data)
        return dataset


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.write(data, 1)
print(memfile.read())

Output:

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

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