I have a shapefile of points that has a field of yield data from harvesting. I am trying to 'smooth' the data by taking each point and collecting all points within 50ft and then update the shapefile with the average of the 50ft. circle. The end goal is that I want a shapefile with every point having a field that contains the average yield of all the other points within 50ft.
Unfortunately python isn't my main programming language and I'm just learning about gis.
The code below is working, but I'm not satisfied with it. What is a performant way to find all points within a given radius of another point, and then do that for every point in the file?
outliers_r is a field on the shapefile where I cleaned up some of the outliers in the data.
yield_smoo is the field I want updated with the average yield.
from osgeo import ogr import fiona import pyproj as proj from shapely.geometry import shape from shapely import geometry #open the file source = ogr.Open("2017_wheat_yield_data.shp", update=True) layer = source.GetLayer() fc = fiona.open("2017_wheat_yield_data.shp") fc_2 = fiona.open("2017_wheat_yield_data.shp") crs_wgs = proj.Proj(init='epsg:4326') crs_bng = proj.Proj(init='epsg:3420') for feat in fc: coords = feat["geometry"]["coordinates"] x, y = coords x1, y1 = proj.transform(crs_wgs, crs_bng, x, y) point_1 = geometry.Point(x1, y1) yield_data =  for feat_2 in fc_2: coords_2 = feat_2["geometry"]["coordinates"] x_2, y_2 = coords_2 x2, y2 = proj.transform(crs_wgs, crs_bng, x_2, y_2) point_2 = geometry.Point(x2, y2) distance = 50 # circle_buffer = point_1.buffer(distance) # if circle_buffer.contains(point_2): # print('circle buffer contains point 2') if point_1.distance(point_2) < distance: data = feat_2["properties"]["outliers_r"] yield_data.append(data) # sum the yields average_yield = sum(yield_data) / float(len(yield_data)) id = int(feat["id"]) fture = layer.GetFeature(id) fture.SetField("yield_smoo", average_yield) layer.SetFeature(fture)