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
up vote 1 down vote accepted

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]
  • I ran your code and while it worked it does not return the Central Feature but the feature closest to the Mean Center. I am looking for a faster way to get the feature with the shortest summed distance to all other features. I also need to spend time getting intimate with numpy :) Your code will come in handy for a rapid way to calculate the Mean Center though which forms another part of my research. Thanks for the input – Clubdebambos Mar 3 '17 at 17:04
  • Absolutely smashing stuff. 3.65 seconds. Exactly what I was looking for. I will break this down using your comments in the code and learn Numpy as I do it. Cheers – Clubdebambos Mar 3 '17 at 17:41
  • I added a change to use numpy's broadcasting, which avoids the nested for loop I had earlier. Always forgetting how to use newaxis... – Logan Byers Mar 3 '17 at 17:43

Your Answer

 
discard

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

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