Algorithm for calculating the Central Feature using open source Python

Question

I want to use an open source programmatic approach through Python for calculating the central feature of a feature class with over 700 points. At the moment the process in the code below works but just very slow. It is just over 2 minutes in comparison to 2 seconds when performed in ArcGIS. The output is the same which is good at least.

The question is; are there any data structures or better algorithms that you can suggest to improve the speed of finding the central feature? I have searched on here and through Google but not returning much information. Has anyone attempted this before? or can anyone put forward an efficient way to find the central feature that trumps what I currently have?

The central feature by definition is the feature whose summed distances to all other features is the shortest.

from osgeo import ogr
from shapely.geometry import Point
from datetime import datetime

start_time = datetime.now()

## set the driver for the data
driver = ogr.GetDriverByName("FileGDB")
## path to the FileGDB
gdb = r"C:\Users\*****\Documents\my_geodata.gdb"
## ope the GDB in write mode (1)
ds = driver.Open(gdb, 1)
## reference the layer using the layers name
lyr = ds.GetLayerByName("my_points")

## massive number that the shortest summed distance sill be less than
shortest__total_distance = 1000000000000.00

## keep track of feature x and y with shortest summed distances
central_x = 0.00
central_y = 0.00

## keep track of the features that have been processed
feature_index = 1

## for each point in the layer
for pnt_from in lyr:
    ## set the total distnace to 0.0
    pnt_total_dist = 0.0
    ## reference the feature at the currnt index
    feature = lyr.GetFeature(feature_index)
    ## access the geometry of that feature
    feature_geom = feature.geometry()
    ## get the x and y coords
    x = feature_geom.GetX()
    y = feature_geom.GetY()
    ## reset the index pointer to the first feature
    lyr.ResetReading()

    ## for each point in the layer
    for pnt_to in lyr:
        ## access the geometry and get the x, y coords
        pt = pnt_to.geometry()
        to_x = pt.GetX()
        to_y = pt.GetY()

        ## calculate the distance between each point and the feature_index point
        pnt_distance = Point(x, y).distance(Point(to_x, to_y))

        ## sum the distances cumulatively for each pair
        pnt_total_dist += pnt_distance

        ## if the total distance goes over the shortest total
        ## then stop for this particular point
        if pnt_total_dist > shortest__total_distance:
            break

    ## if the total distance is shorter than the current shortest total
    ## it becomes the new shortest total and updates the central x and y
    if pnt_total_dist < shortest__total_distance:
        shortest__total_distance = pnt_total_dist
        central_x = x
        central_y = y

    ## reset the index pointer to the first feature
    lyr.ResetReading()
    ## increase the feature index
    feature_index += 1
    print feature_index
    ## when the feature index is > the feature count the process is over
    if feature_index > lyr.GetFeatureCount():
        break

## print the coordinates of the central feature
print central_x, central_y

## This part creates a feature class with the central feature as a point.
"""
## create a new point layer with the same spatial ref as lyr
out_lyr = ds.CreateLayer("central_feature", lyr.GetSpatialRef(), ogr.wkbPoint)

## define and create new fields
x_fld = ogr.FieldDefn("X", ogr.OFTReal)
y_fld = ogr.FieldDefn("Y", ogr.OFTReal)
out_lyr.CreateField(x_fld)
out_lyr.CreateField(y_fld)

## create a new point for the mean center
pnt = ogr.Geometry(ogr.wkbPoint)
pnt.AddPoint(central_x, central_y)

## add the central feature to the new layer
feat_dfn = out_lyr.GetLayerDefn()
feat = ogr.Feature(feat_dfn)
feat.SetGeometry(pnt)
feat.SetField("X", central_x)
feat.SetField("Y", central_y)
out_lyr.CreateFeature(feat)

feat = None
"""

ds = None

print datetime.now() - start_time

Answer 1

EDIT:

Create a NxN array of distances and sum them along an axis. Get the min summed distance.

from osgeo import ogr
import numpy as np      # USING NUMPY 

## set the driver for the data
driver = ogr.GetDriverByName("FileGDB")
## path to the FileGDB
gdb = r"C:\Users\*****\Documents\my_geodata.gdb"
## open the GDB in write mode (1)
ds = driver.Open(gdb, 1)
## reference the layer using the layers name
lyr = ds.GetLayerByName("my_points")

# get numpy array of (x,y)'s for every point
xy_arr = np.ndarray((len(lyr), 2), dtype=np.float)
for i, pt in enumerate(lyr):
    g = pt.geometry()
    xy_arr[i] = (g.GetX(), g.GetY())

# construct NxN array of inter-point distances
# numpy broadcasting is used here
pt_dist_arr = np.linalg.norm(a-a[:,np.newaxis, :], axis=2)

# sum distances for each point
summed_distances = np.sum(pt_dist_arr, axis=0)

# index of point with minimum summed distances
index_central_feat = np.argmin(summed_distances)

# position of the point with min distance
central_x, central_y = xy_arr[index_central_feat]

about

I guess I was assuming istropic spatial variance in my original post.

original:

Your algorithm is worst case O(n^2) I think, but the real hits might be the lookup time for the points within the OGR structure and the many data type conversions. lyr is being both iterated over (the pointer to the current feature is being moved serially) and also random access is occuring (lyr.GetFeature()). All while making a bunch of one-time objects.

My immediate thought for this task was to find the point nearest to the centroid of all the points. Should be O(n) or better.

from osgeo import ogr
import numpy as np      # USING NUMPY 

## set the driver for the data
driver = ogr.GetDriverByName("FileGDB")
## path to the FileGDB
gdb = r"C:\Users\*****\Documents\my_geodata.gdb"
## ope the GDB in write mode (1)
ds = driver.Open(gdb, 1)
## reference the layer using the layers name
lyr = ds.GetLayerByName("my_points")

# get numpy array of (x,y)'s for every point
xy_arr = np.ndarray((len(lyr), 2), dtype=np.float)
for i, pt in enumerate(lyr):
    g = pt.geometry()
    xy_arr[i] = (g.GetX(), g.GetY())

# get average x,y (centroid of all points)
avg_x, avg_y = np.mean(xy_arr, axis=0)
print "centroid", avg_x, avg_y

# 1-dimensional array of distances of point to average
distances = np.linalg.norm(xy_arr - [avg_x, avg_y], axis=1)

# index of point with minimum distance to centroid
index_central_feat = np.argmin(distances)

# position of the point with min distance
central_x, central_y = xy_arr[index_central_feat]