CODEWITHSUNDEEP
pythonrgeopandasterrarasterio
89_Simple
89_Simple

Reputation: 3805

area weighted average of raster using polygon - R vs Python

I have an R function in terra package that takes a raster and polygon and calculates the average of the raster for the polygon. In doing so, it uses the area intersection of each raster cell with the polygon as a weight to do the averaging (for e.g. some raster on polygon boundary could be partially intersecting with the polygon). The weights = T argument does this for me. Here's the R code:

library(terra)
terra::rast(my_raster, my_shp, fun = 'mean', na.rm = T, weights = T, touches = T)

I want to do the equivalent in python:

import rasterio
import geopandas as gpd
import numpy as np 

out_image, out_transform = rasterio.mask.mask(my_raster, my_shp, crop=True, all_touched=True)
mean_val = np.mean(out_image) # Calculate the mean value of the masked raster

However, I am not able to find an argument in python that accounts for raster cell's intersection area with the polygon to calculate the average. How do people go about doing this in python?

Upvotes: 1

Views: 472

Answers (1)

Sébastien Langlois
Sébastien Langlois

Reputation: 21

If you work with the xarray format in python, then you can do it with the clisops library. Specifically, lookup for create_weight_masks function.

import geopandas as gpd  
import xarray as xr  
from clisops.core.subset import create_weight_masks  
ds = xr.open_dataset(path_to_tasmin_file)  
polys = gpd.read_file(path_to_multi_shape_file)  

# Get a weight mask for each polygon in the shape file
mask = create_weight_masks(ds_in=ds, poly=polys)

Upvotes: 2

Related Questions