1

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?

  • 1
  • I can see WriteArray (to be used with gdal.band type of data) does what I'm looking for. Any idea whether there's anything similar for raster data opened with 'rasterio'? The reason I'm asking is because the processing done so far is with rasterio (like the use of 'mask' for example). If I open the raster with gdal instead I guess I would need to find a 'gdal' way to do the mask... – Pitrako Junior Jul 18 at 14:42
  • I haven't used rasterio but as far as I know, it runs on top of gdal. That means that there is a chance their .write() method accepts an offset as well. You could take a look at the documentation or tweak around with the .write() method to see if its possible. – Marcelo Villa Jul 18 at 14:53
0

I would recommend using rioxarray.

import numpy
import rioxarray
from shapely.geometry import mapping

xds = rioxarray.open_rasterio("NI_RIVER_ALR_1000.tif").sel(band=1).squeeze()
gdb_vml_polygs = 'RSA_VML_BUILDINGS_NI.gdb'
fc_vml_polygs = 'VML_Buildings_NI'
vml_polygs = gpd.read_file(gdb_vml_polygs, driver='FileGDB', layer=fc_vml_polygs)
out_xds = xds.copy()
for geom in vml_polygs.geometry.apply(mapping):
    clipped = xds.rio.clip([geom], vml_polygs.crs, drop=False)
    out_xds = out_xds.where(numpy.isnan(clipped), clipped.mean())
out_xds.rio.to_raster("out_averages.tif")

I tested it on this input raster:

enter image description here

And with four regions to average over, this was the result:

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.