4

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)
  • Welcome to GIS SE. Please pare this down to a single, focused question. Once that is complete, vote to reopen. – Aaron Aug 8 '17 at 4:57
  • @Aaron I know that the question is a bit broad. My knowledge is weak enough in this area that I wanted to give enough information that I could be pointed in the right direction. xunilk's entire answer is very helpful, but especially the statement "I think that fiona was a good choice..." was the type of answer I wanted. At this point, my only question is how to write the calculated data back to the same file with fiona. Do you have advice on how to reword the question without making xunilk's answer irrelevant? – javier Aug 8 '17 at 14:31
2

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.

enter image description here

  • Can fiona update an existing shapefile or do I need to write the calculated data to a new shapefile? – javier Aug 8 '17 at 14:32
  • I think that it's preferable to write the calculated data to a new shapefile. See this answer: gis.stackexchange.com/questions/215963/… – xunilk Aug 8 '17 at 14:46
  • Please, edit this question to fit one topic. I would vote to reopen. – xunilk Aug 8 '17 at 14:50
  • Edited it. Did I make it clear enough now? – javier Aug 8 '17 at 15:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.