5

I have a raster with a width of 2979 and a height of 1867 pixels. I need to split the raster into tiles of equal size (e.g. 5 kilometres x 5 kilometres). To do so, I need to update the height and width of the raster to fit an exact (i.e. integer) number of tiles. If the chosen dimension is 5km x 5km and the spatial resolution of the raster is 10 metres, the width of the raster would need to be increased to 3000 pixels and the height to 2000 pixels. This would create 6 (3000 / 500) x 4 (3000 / 500) = 24 tiles. I am using rasterio to work with the raster. This is my current code:

# Determine padding
padding_x = int(new_raster_width - raster.width)
padding_y = int(new_raster_height - raster.height)

# Update values based on padding
out_meta.update({
    'width': new_raster_width,
    'height': new_raster_height
})

# Loop through raster bands
for band_no in range(1, raster.count + 1):
    # Read band
    raster_band = raster.read(band_no)

    # Add padding
    padded_raster_band = np.pad(
        raster_band, 
        pad_width = (padding_x, padding_y), 
        mode = 'constant',
        constant_values = 0)

    # Append padded raster band to list
    padded_bands.append(padded_raster_band)

with rio.open(raster_out, 'w', **out_meta) as dest:
    for band_nr, src in enumerate(padded_bands, start=1):
        dest.write(src, band_nr)

I am trying to do the following: Calculate by how much pixels I need to increase the raster, update the width and height, loop through each band of the raster and update the values and, finally, write these bands into a new dataset.

This works. I do have, however, the following problem: The raster is also scaled, which results in an incorrect 'projection'. I have also included two images of the raster before and after increasing the size (although it might be hard to tell from the images).

How do I add a padding/increase the raster size without changing the scale of the raster?

enter image description here

enter image description here

1
  • Certainly you do not need to edit the source file that way just for cutting tiles but you can define tile extents which are partly or even totally outside the extent of the source image. However, I do not know how to do that with rasterio. With gdal_translate it can be done with -srcwin gdal.org/programs/…. For example gdal_translate -srcwin -100 -150 200 200 input.tif output.tif writes a 200x200 sixed tile that has the top left corner 100 pixels left and 150 pixels upwards from the top left corner of the source image.
    – user30184
    Oct 29, 2021 at 12:38

2 Answers 2

4

You don't need to pad at all. Just read each tile, using boundless=True and let rasterio pad for you.

from itertools import product
import rasterio as rio
from rasterio.windows import Window

raster = '/path/to/input.tif'
width, height = 500, 500  # size of tiles in pixels
fill_value = 0

with rio.open(raster) as src:
    offsets = product(range(0, src.meta['width'], width), range(0, src.meta['height'], height))
    for col_off, row_off in offsets:
        window = Window(col_off=col_off, row_off=row_off, width=width, height=height)
        data = src.read(boundless=True, window=window, fill_value=fill_value)
        print(data.shape)
3

If you change the origin of the raster (typically the upper left) by adding the padded values, you need to also change the affine transformation (src.transform). But if you only add padding to the bottom and right of your original extent, everything should theoretically stay where it should be.

However as far as I can tell you don't exactly pad your array the right way. The way you have it right now I think it pads both axis the same amount. Try it this way and see if it works:

# Add padding
padded_raster_band = np.pad(
    raster_band, 
    pad_width = ((0, padding_y), (0, padding_x)), 
    mode = 'constant',
    constant_values = 0)

The first tubple of pad_width is the padding of the y axis, with the first value being the padding before the array and second after the array. The same is the case for the second tuple only for the x axis.

0

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.