Python collect all points that meet a criteria for every point in the shapefile

Question

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?

The outliers_r is a field on the shapefile where I cleaned up some of the outliers in the data.

The 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)

Answer 1

You have many questions here so, I'm only going to answer question related with fiona and distances. I think that fiona was a good choice to open your file and you don't need to open it three times. On the other hand, for avoiding measure distances twice, you can use itertools python module.

To try out my approach, I generated 100 arbitrary points in Kansas South (CRS 3420) and all distances (into list comprehensions) were measured with next code:

import fiona
from shapely.geometry import Point
import itertools

fc = fiona.open("/home/zeito/pyqgis_data/2017_wheat_yield_data.shp")

points = [ Point(feat["geometry"]["coordinates"]) for feat in fc ]

comb = range(len(points))

distances = [ [i, j, points[i].distance(points[j])] 
              for i, j in itertools.combinations(comb, 2) ]

distances_less50 = [ [i, j, points[i].distance(points[j])]
                     for i, j in itertools.combinations(comb, 2) 
                     if points[i].distance(points[j]) <= 50 ]

Then, in each case, first list distance is:

>>>distances[0]
[0, 1, 1344.1704643424168]
>>>distances_less50[0]
[5, 54, 35.45937075954941]
>>>len(distances)
4950
>>>len(distances_less50)
6 

where it can be observed its respective indexes.

At next image, it can be visualized first pair where distance is lesser than 50 feet.