Skip to main content
4 of 8
Title edit.
TomazicM
  • 27k
  • 23
  • 31
  • 42

Updating a subset of raster values using open source Python (rasterio, for instance)?

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?

Pitrako Junior
  • 1.3k
  • 18
  • 31