# Performance problem with getting average pixel values within buffered circles

I have a EPSG:4326 raster file with pixel size of 15 arcseconds, and a list of longitudes and latitudes. For each pair of coordinate, I want to create a buffered circle of X kilometers surrounding them, and aggregate over the pixel values within each circle.

Following this thread: Rasterio. Pixel value of a point X km away from current point:

1. I have used Shapely to get a circle of defined radius from a centre point (after reprojecting the EPSG:4326 file to a projected coordinate system).
2. Then, I use `rasterio.features.geometry_mask` to calculate a numpy mask given the Shapely circle.
3. Finally, I calculate the average pixel value on the masked array with `np.ma.mean`

The problem is the `np.ma.mean` function is too slow for my numpy array (it has the shape of `(34779, 45217)`, and dtype of `float32`), and it can take hours to perform the operation over my list of 5000 longtitudes and latitudes. Is there any other way? I am open to any tool, including Python, R, and any and all software.

Below is the annotated code for your reference.

``````from functools import partial
import pyproj
from shapely.geometry import Point
from shapely.ops import transform
from rasterio.features import geometry_mask
import rasterio
import pandas as pd
import numpy as np

sample = pd.read_csv(lat_lon_url) # read in dataframe that contains longtitude and latitude
raster = rasterio.open(raster_url)

'''
apply the get_mean function for each row of the sample dataframe.
each row of the sample dataframe has a longitude and latitude.
'''
sample['lightscore'] = sample.apply(lambda row: get_mean(row['longtitude'], row['latitude'],
raster.height, raster.width, raster.transform,
band, radius = 5000), axis = 1)

``````

Below is the `get_mean` function.

``````'''
get_mean function
1. creates a circle with specified radius surrounding the point with longitude & latitude
2. aggregate over the circle to get average pixel value
'''

def get_mean(lon, lat,
raster_height, raster_width, raster_transform,
band, radius: float) -> float:

'''
Define the projected and geographic coordinate system
'''
local_azimuthal_projection = "+proj=aeqd +R=6371000 +units=m +lat_0={} +lon_0={}".format(
lat, lon
)
wgs84 = "+proj=longlat +datum=WGS84 +no_defs"

'''
Define the projection.
'''
wgs84_to_aeqd = partial(
pyproj.transform,
pyproj.Proj(wgs84),
pyproj.Proj(local_azimuthal_projection),
)
aeqd_to_wgs84 = partial(
pyproj.transform,
pyproj.Proj(local_azimuthal_projection),
pyproj.Proj(wgs84),
)

center = Point(float(lon), float(lat))
point_transformed = transform(wgs84_to_aeqd, center) # reproject the point from geographic to projected system
circle_poly = transform(aeqd_to_wgs84, buffer) # reproject the buffered circle from projected back to geographic system

'''
Calculate a numpy mask given the Shapely circle.
'''
out_shape = (raster_height, raster_width),
transform = raster_transform) # this takes a while for each record of my sample dataframe,
# but not as long as the np.ma.mean below

'''
Create masked numpy array and aggregate over it
'''
mean = masked_data.mean() # this takes EXTREMELY LONG for each record of my sample dataframe.

return mean
``````
• Intuitively, I think that running `mean` on a `np.ma.array` should be quite fast. What is the shape and dtype of your array? Commented Jan 20, 2022 at 9:33
• @StefanBrand I have updated my question to include the shape and dtype of my numpy array. Its shape is `(34779, 45217)`, and the dtype is `float32` Commented Jan 20, 2022 at 14:41
• If I'm not mistaken, your array has a size of ~6 GiB: `GiB = 34779 * 45217 * 32 / 8 / 1024 / 1024 / 1024`. Trying to allocate a random numpy array of this shape crashes here. Commented Jan 20, 2022 at 15:01
• I'd probably generate a geodata file for all circles and then use `rasterstats` or even QGIS to do the zonal statistics for me. They are surely better suited for such large datasets than looping with numpy. Commented Jan 20, 2022 at 15:16
• @StefanBrand I have not used these tools, but from your text, I assume that QGIS is more suitable for larger datasets compared to `rasterstats`, and both of them are definitely better than `numpy`? Commented Jan 20, 2022 at 16:18

PyQGIS manages to perform this task significantly faster than `numpy` or `rasterstats`. Here is the link to the cookbook: https://docs.qgis.org/3.16/en/docs/pyqgis_developer_cookbook/index.html.

To complete the task:

... a EPSG:4326 raster file with pixel size of 15 arcseconds, and a list of longitudes and latitudes. For each pair of coordinate, ... create a buffered circle of X kilometers surrounding them, and aggregate over the pixel values within each circle.

You can use the below template:

``````from qgis.core import *
from qgis.PyQt.QtCore import QVariant
from qgis.analysis import QgsZonalStatistics

import pandas as pd

# Inits app
QgsApplication.setPrefixPath(path_to_qgis, True)
qgs = QgsApplication([], False)
qgs.initQgis()

# Sets up original CRS
epsg4326 = QgsCoordinateReferenceSystem("EPSG:4326")
transformContext = QgsProject.instance().transformContext()

# Creates shapefile to write buffer into
filename = "buffers.shp"

fields = QgsFields()
fields.append(QgsField(lon_col_name, QVariant.Double))
fields.append(QgsField(lat_col_name, QVariant.Double))

save_options = QgsVectorFileWriter.SaveVectorOptions()
save_options.driverName = "ESRI Shapefile"
save_options.fileEncoding = "UTF-8"

writer = QgsVectorFileWriter.create(
filename,
fields,
QgsWkbTypes.Polygon,
epsg4326,
transformContext,
save_options
)

# Loads .csv file as vector layer
vlayer = QgsVectorLayer(data_url, "coordinates", "ogr")
features = vlayer.getFeatures()

field_names = [ field.name() for field in vlayer.fields() ]
lat_col_idx = field_names.index("lat")
lon_col_idx = field_names.index("lon")

for f in features: # each feature is a point
attributes = f.attributes()
lat = float(attributes[lat_col_idx])
lon = float(attributes[lon_col_idx])

# Sets up projected CRS and transformer
local_azimuthal = QgsCoordinateReferenceSystem()
proj4_str = "+proj=aeqd +R=6371000 +units=m +lat_0={} +lon_0={}".format(lat, lon)
local_azimuthal.createFromProj(proj4_str)

wgs84_to_azimuthal = QgsCoordinateTransform(epsg4326, local_azimuthal, transformContext)
azimuthal_to_wgs84 = QgsCoordinateTransform(local_azimuthal, epsg4326, transformContext)

# Projects the point, creates buffer, and reprojects to original CRS
projected_point = wgs84_to_azimuthal.transform(QgsPointXY(lon, lat))
point_geom = QgsGeometry.fromPointXY(projected_point)
radius = 5000 # in meters
segment = 18
buffer = point_geom.buffer(radius, segment)
buffer.transform(azimuthal_to_wgs84)

# Writes buffer to file
fet = QgsFeature()
fet.setGeometry(buffer)
fet.setAttributes([lon, lat])

# flush file to avoid memory leak
del writer

# Reads in the buffer shape files
filename = "buffers.shp"
vlayer = QgsVectorLayer(filename, "buffers", "ogr")

# Calculates zonal mean-s for each raster
raster = QgsRasterLayer(raster_url)
zoneStats = QgsZonalStatistics(vlayer, raster, stats = QgsZonalStatistics.Statistics(QgsZonalStatistics.Mean))
zoneStats.calculateStatistics(None)

# Writes zonal mean-s to csv