I have a source raster data set (tif file) and a number of polygon features. For each one of those polygons I need to: a) 'clip' the raster. b) change values of those clipped raster cells (with the average value). c) write new values back to the source raster. I've managed to do a and b, which means I now have a `numpy.ndarray` variable where all cells have been populated with the average value. How can I write this back to the source raster? I believe I'd be able to write a script which could do such a thing. I have arrays dimensions (width and height) and top left cell's position, so I guess it'd just be a matter of operating with these arrays. However, I was wondering whether there's any existing function that does that already. My code: import os, sys, datetime, time import geopandas as gpd import gdal import rasterio from rasterio.mask import mask from fiona.crs import from_epsg import numpy as np import pycrs alr_path = r'Z:\GRAU_Team_Admin\Alfonso Jimenez\Flood Modelling\Building Footprint\OUTPUT\NI_RIVER_ALR_1000.tif' gdb_vml_polygs = r'\\lwukwvdi11\data\Data\UK\UK_Mapping\VML\data\Europa20190403\RSA_VML_BUILDINGS_NI.gdb' fc_vml_polygs = r'VML_Buildings_NI' alr = rasterio.open(alr_path) vml_polygs = gpd.read_file(gdb_vml_polygs, driver='FileGDB', layer=fc_vml_polygs) for index, row in vml_polygs.iterrows(): #row[3] contains the geometry window, out_transform = mask(alr, row[3], all_touched=True, crop=True) if np.all([window < 0]): avg = -1 else: avg = window[window != -1].mean().item() window_avg = (np.where(window!=-1, int(round(avg,0)), window)) In the script `window_avg` is a numpy array containing the new values. How can I write these values back into the source raster?