0

I am trying to create a large raster of 64-bit integers from the metadata of another file. This is done row by row, writing the output with the rasterio Window class. Below is a very simplified example that assigns the column number to each cell:

import rasterio
from rasterio.windows import Window
import numpy as np

input_path = "/path/to/my/large/raster.vrt"
output_path = "/path/to/the/output.gtiff"

input = rasterio.open(input_path)

out_file = rasterio.open(
        output_path,
        'w',
        driver='GTiff',
        height=input.height,
        width=input.width,
        count=1,
        dtype=np.int64,
        crs=input.crs,
        transform=input.transform
    )

for j in range(input.height):

    array = np.full((input.width, 1), 10000).astype(np.int64)

    for i in range(input.width):
        array[i, 0] = i

    window = Window(0, j, array.shape[0], array.shape[1])
    out_file.write(array, window=window, indexes=1)

out_file.close()

This runs fine, but the output is a zeroed raster:

gdalinfo -mm /path/to/the/output.gtiff
[...]
Band 1 Block=89057x1 Type=Int64, ColorInterp=Gray
    Computed Min/Max=0.000,0.000

Looks like no array is ever written to the file, or that something happens to the integers sent for writing. Any clues on what I may be missing?

1
  • @user2856 could you add that as an answer? Commented Mar 14 at 14:41

1 Answer 1

1

You have rows and cols reversed, use

array = np.full((1, input.width), 10000).astype(np.int64) 

and

array[0, i] = i

But note that looping through every pixel individually is a really inefficient way of doing things.

For your trivial writing column number example, you could use np.arange by rows only:

array = np.arange(input.width, dtype=np.int64).reshape(1, input.width)
for j in range(input.height):
    window = Window(0, j, input.width, 1)
    out_file.write(array, window=window, indexes=1)

Or in one go using array broadcasting

cols = np.arange(input.width, dtype=np.int64).reshape(1, input.width)
rows = np.zeros((input.height, 1), dtype=np.int64)
array = rows + cols
out_file.write(array, indexes=1)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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