I am trying to read a raster with overlapping windows and doing some calculations on those windows. I am using the code from this answer. I want to do calucations on overlapping windows but want to write the windows without overlap.
from itertools import product
import rasterio
from rasterio import windows
def overlapping_windows(src, overlap, width, height, boundless=False):
""""width & height not including overlap i.e requesting a 256x256 window with
1px overlap will return a 258x258 window (for non edge windows)"""
offsets = product(range(0, src.meta['width'], width), range(0, src.meta['height'], height))
big_window = windows.Window(col_off=0, row_off=0, width=src.meta['width'], height=src.meta['height'])
for col_off, row_off in offsets:
window = windows.Window(
col_off=col_off - overlap,
row_off=row_off - overlap,
width=width + overlap * 2,
height=height + overlap * 2)
if boundless:
yield window, transform
else:
yield window.intersection(big_window), transform
def process_image(src_img, dst_img, band_id=1):
with rasterio.open(src_img) as src:
kwargs = src.meta
with rasterio.open(dst_img, 'w', **kwargs) as dst:
for window, transform in overlapping_windows(src, overlap, width, height, boundless=False):
src_data = src.read(band_id, window=window)
dst_data = src_data ** 2 # Do the Processing Here
# write raster
profile = {
'driver': 'GTiff',
'tiled': True,
'compress': None
}
dtype = rasterio.dtypes.get_minimum_dtype(dst_data)
crs = CRS.from_epsg(epsg)
height = dst_data.shape[1]
width = dst_data.shape[2]
count = dst_data.shape[0]
with rasterio.Env():
profile.update(
dtype=dtype,
count=count,
height=height,
width=width,
crs=crs,
transform=transform
)
with rasterio.open(out_path, 'w', **profile) as dst:
for band_idx in range(count):
arr_out = dst_data[band_idx, :, :, 0].astype(rasterio.dtypes.get_minimum_dtype(dst_data))
dst.write(arr_out, 1 + band_idx)
return 0
This writes the overlapping windows, while I want to write non-overlapping windows but could not get my head around how to that.
