I am using this questions solution to split a raster I have into 1500x1500 tiles. The code works. I have recently been trying to figure out how to set the starting and ending points of the raster based on the extent of a shapefile. For example, the cat picture attached, the red box in the shapefile extent I want to use as the window for tiling. The blue boxes are the 1500x1500 tiled images. The white part makes up the additional pixels that are outside of the shapefiles extent but still reach the 1500x1500 size.
I tried changing the image's meta data 'transform' to the shapefiles extent through the commented out code below, but get TypeError: 'Affine' object does not support item assignment. So it seems that I cannot change it that way. I thought about clipping the image to the shapefile and then running with that. But I'd like to bypass that step since clipping large images (10GB) can take a while.
Any suggestions?
import os
from itertools import product
import rasterio as rio
from rasterio import windows
import fiona
from shapely.geometry import shape
out_path = 'C:/Users/z/Desktop/output'
output_filename = 'item'
Set the wdith and height of wanted split image size
width = 1500
height = 1500
Function to split the raster
def get_tiles(ds, width=width, height=width):
nols, nrows = ds.meta['width'], ds.meta['height']
offsets = product(range(0, nols, width), range(0, nrows, height))
big_window = windows.Window(col_off=0, row_off=0, width=nols, height=nrows)
for col_off, row_off in offsets:
window = windows.Window(col_off=col_off, row_off=row_off, width=width, height=height).intersection(big_window)
transform = windows.transform(window, ds.transform)
yield window, transform
training = fiona.open(shapefile.shp)
x = next(iter(training))
extent = shape(x['geometry']).bounds
xmin = extent[0]
ymin = extent[1]
i = 0
with rio.open('C:/Users/z/Desktop/ortho/Orthomosaic.tif') as inds:
# inds.meta['transform'][2] = xmin
# inds.meta['transform'][5] = ymin
tile_width, tile_height = width, height
meta = inds.meta.copy()
for window, transform in get_tiles(inds):
meta['transform'] = transform
meta['width'], meta['height'] = window.width, window.height
outpath = os.path.join(out_path, output_filename + str(i) + '.tif')
i = i + 1
# with rio.open(outpath, 'w', **meta, compress='DEFLATE') as outds:
with rio.open(outpath, 'w', **meta) as outds:
outds.write(inds.read(window=window))
