How to rasterize point data using average value per cell (not replace or add)

[dluks]

Hey all, I have a geopandas.GeoDataFrame containing point observations that I would like to rasterize into a grid with a certain resolution, such that, for cells that contain multiple observations, the resulting grid cell value is the average value of the contained points.

In an ideal world I'd use geocube.api.core.make_geocube due to its ease of use. My current code looks something like this.

import json

import geopandas as gpd
import numpy as np
from geocube.api.core import make_geocube
from shapely.geometry import box, mapping, Point

# global bounding box
bbox = (-180, -90, 180, 90)

# some random points
points = [Point(np.random.uniform(-180, 180), np.random.uniform(-90, 90)) for _ in range(1000)]

# some random data values
data = np.random.rand(1000)

# create GeoDataFrame
gdf = gpd.GeoDataFrame({'data': data}, geometry=points, crs='EPSG:4326')

grid = make_geocube(
    vector_data=gdf,
    measurements=[column],
    resolution=(-0.5, 0.5),
    geom=json.dumps(mapping(box(*bbox))),
) 

The problem here, if I understand correctly, is that geocube's default rasterize_function (geocube.rasterize.rasterize_image) utilizes rasterio.enums.MergeAlg as the only options for merging point values that are contained in a single grid cell, i.e. GDAL's REPLACE or ADD methods. Additionally, in my case I don't want to do any interpolation of missing values, i.e. the resulting grid should be sparse (I realize my example code doesn't really simulate that part).

Is there a way to easily implement a custom AVERAGE merge algorithm in this case, or will I need to do it all manually, e.g. with rioxarray or other tools? And in case I'm missing something, is there another built-in method geocube that can do this while avoiding interpolation of the entire grid?

[snowman2]

I am not aware of any. However, as a workaround, you could use ADD and rasterize two vectors. One with the values and a duplicate with a 1 replacing the values. Then, divide the raster with values by the raster generated by the vector filled with ones.

[dluks]

OK, that could work, though it feels a bit circuitous. I'm a bit surprised at this limitation (in this case I guess the limitation goes all the way back to GDAL) since this is a pretty common operation. For example, taking a set of point observations such as land surface temperature and simply gridding them without interpolating across the entire grid extent. ADD wouldn't make any sense in that case, and REPLACE is just throwing away data.

Any thoughts on how this might be incorporated into geocube and would a PR be welcome? Either into geocube.rasterize.rasterize_image or geocube.rasterize.rasterize_points_griddata?

[snowman2]

Any thoughts on how this might be incorporated into geocube

• The best case scenario would be to add it to GDAL
• The second best case scenario would be to add it to rasterio
• As a last resort, it could be added to geocube. You could accomplish this in rasterize_image.