1

I have lots of start points and a handful of destination points. Using 'Point Distance' from the toolbox takes ages to calculate and write the results. Hence, I was looking for alternatives and went the Numpy, SciPy way.

Here how I calculate the distances, and it only takes some seconds if at all.:

import arcpy
import scipy.spatial
arcpy.env.workspace = 'Drive:/Path/To/Database.gdb'
start = arcpy.da.FeatureClassToNumPyArray('Foo', ['Bar', 'SHAPE@XY'])
dest = arcpy.da.FeatureClassToNumPyArray('John', ['Doe', 'SHAPE@XY'])
arr_dist = scipy.spatial.distance.cdist(start['SHAPE@XY'], dest['SHAPE@XY'], 'euclidean')

That gives me in my case a matrix with the dimension of (670911, 39). That means, each row is representing a feature from start, while each column holds the distance for the particular point from the start FeatureClass to all destination points. The order of the destination follows the order of OIDs.

My question is: How do I get the information from 'Bar' and 'Doe' to the arr_dist and get it exported into a table in my database.

The arcpy.da.NumPyArrayToTable function requires a structured array and the dtype for every column, which seems to be pretty cumbersome considering the number of columns I have.

[EDIT]

I tried to append the information from 'Doe' to arr_dist using the following line:

import numpy.lib.recfunctions
arr_doe = numpy.lib.recfunctions.append_fields(arr_dist.transpose(), 'doe', dest['Doe'])

That has the reproducible effect that ArcGIS (10.4) crashes. When running it in IDLE a 'MemoryError' is returned.

0

Here is an example and it uses rec_append_fields from recfunctions. Most of the code is a sample data preparation section and formatting. I included it since you might be interested in that as well. Generating the names for the output fields is fairly straightforward should you not have a large number to generate.

"""
Script:   append_recarray_cdist_example.py
Author:  Dan.Patterson@carleton.ca
"""
from scipy.spatial.distance import cdist
import numpy as np
from numpy.lib import recfunctions as rfn
formatter = {'float': '{:0.3f}'.format,
             'float64': '{:0.3f}'.format }
np.set_printoptions(edgeitems=2,linewidth=80,precision=2,suppress=True,
                    threshold=6, formatter=formatter)
#
# ---- data example and preparation
Xs = [0,2,1,3,5,4,7,6,9,8]
Ys = [9,7,6,3,4,1,5,8,2,0]
N = len(Xs)
IDs = [i for i in range(N)]
foo_data = list(zip(IDs,zip(Xs,Ys)))
bar_data = [(IDs[i],[Xs[i]+10,Ys[i]+10]) for i in range(N)]
dt1 = [('Foo', '<i4'), ('XY', ('<f8', (2,)))]
dt2= [('Bar', '<i4'), ('XY', ('<f8', (2,)))]
orig = np.array(foo_data,dtype=dt1)
dest = np.array(bar_data, dtype=dt2)
# --- now some work
dist_arr = cdist(orig['XY'], dest['XY'], 'euclidean')
names = ["Bar{}".format(i) for i in range(N)]
final = rfn.rec_append_fields(orig, names, dist_arr) #name corrected
frmt = """
From array...
{}
To array...
{}
Distance array sample...
{!r:}
....snip...
{!r:}
"""
p0 = (final[['Foo','Bar0']]).reshape((N,1))
p1 = (final[['Foo','Bar9']]).reshape((N,1))
print(frmt.format(orig, dest, p0, p1))

The results are as follows

>>> 
From array...
[(0, [0.0, 9.0]) (1, [2.0, 7.0]) ..., (8, [9.0, 2.0]) (9, [8.0, 0.0])]
To array...
[(0, [10.0, 19.0]) (1, [12.0, 17.0]) ..., (8, [19.0, 12.0]) (9, [18.0, 10.0])]
Distance array sample...
array([[(0, 14.142135623730951)],
       [(1, 14.422205101855956)],
       ..., 
       [(8, 19.235384061671343)],
       [(9, 18.027756377319946)]], 
      dtype=[('Foo', '<i4'), ('Bar0', '<f8')])
....snip...
array([[(0, 19.1049731745428)],
       [(1, 17.46424919657298)],
       ..., 
       [(8, 16.278820596099706)],
       [(9, 14.142135623730951)]], 
      dtype=[('Foo', '<i4'), ('Bar9', '<f8')])

Additions: A bit of extra information. Structured arrays and recarrays have no knowledge of columns because they are structured and the columns are accessed by the column name. So in the case of the 'final' array, there are 10 records and each column has a name. To access the information for a particular column or columns you need to employ standard array slicing

>>> final.shape
(10,)
>>> final.dtype.names
('Foo', 'XY', 'Bar0', 'Bar1', 'Bar2', 'Bar3', 'Bar4', 'Bar5', 'Bar6', 'Bar7', 'Bar8', 'Bar9')
>>> 
>>> final['Bar0']
array([14.142, 14.422, ..., 19.235, 18.028])
>>> final[['Foo','Bar0']]
array([(0, 14.142135623730951), (1, 14.422205101855956), ...,
       (8, 19.235384061671343), (9, 18.027756377319946)], 
      dtype=[('Foo', '<i4'), ('Bar0', '<f8')])
>>> 

ADDENDUM 2

A demo of appending two arrays columnwise

# coding: utf-8
"""
Script:  recarray_append_flds.py
Author:  Dan.Patterson@carleton.ca
Purpose: Demonstrate appending columns to a structured/recarray
"""
import numpy as np
from numpy.lib import recfunctions as rfn

# make a reference point file
n = 10
dt0 = [('Id','<i4'),('xy','<f8', (2,))]
xy = np.arange(20).reshape(n,2)
a = np.zeros((n,),dtype=dt0)
a['Id'] = np.arange(n).tolist()
a['xy'] = xy
#
# make some values for combining
#   Now, this is out there, some fancy stuff.  Make 3 float fields (m).
#   The following handles the formatting and empty array construction with default
#   field naming. defaultfmt= "A, B, C" would also be acceptable
m = 3
f = ",".join(['f8' for i in range(m)])  # for the common float fields
dt1 = np.lib._iotools.easy_dtype( f, defaultfmt="A_%02i")
names = list(dt1.names)                 # keep the names for later on
r = np.random.randint(1,5,(n,m))        # if you want integers in n*m format
r = r.astype('float')                   # or as floats in an n*m array
# perform the combination of arrays
b =rfn.append_fields(a,names,np.hsplit(r,3),usemask=False)
frmt = """
Input array
{}
Secondary data
{}
Extended array
{}
"""
print(frmt.format(a.reshape(n,1),r,b.reshape(n,1)))

Results in the following which can be sent back to ArcMap using NumPyArrayToFeatureClass

Input array
[[(0, [0.0, 1.0])]
 [(1, [2.0, 3.0])]
 [(2, [4.0, 5.0])]
 [(3, [6.0, 7.0])]
 [(4, [8.0, 9.0])]
 [(5, [10.0, 11.0])]
 [(6, [12.0, 13.0])]
 [(7, [14.0, 15.0])]
 [(8, [16.0, 17.0])]
 [(9, [18.0, 19.0])]]
Secondary data
[[ 4.  4.  3.]
 [ 1.  4.  2.]
 [ 2.  4.  1.]
 [ 2.  3.  2.]
 [ 4.  3.  4.]
 [ 2.  4.  4.]
 [ 4.  1.  2.]
 [ 4.  4.  2.]
 [ 2.  2.  3.]
 [ 1.  2.  3.]]
Extended array
[[(0, [0.0, 1.0], 4.0, 4.0, 3.0)]
 [(1, [2.0, 3.0], 1.0, 4.0, 2.0)]
 [(2, [4.0, 5.0], 2.0, 4.0, 1.0)]
 [(3, [6.0, 7.0], 2.0, 3.0, 2.0)]
 [(4, [8.0, 9.0], 4.0, 3.0, 4.0)]
 [(5, [10.0, 11.0], 2.0, 4.0, 4.0)]
 [(6, [12.0, 13.0], 4.0, 1.0, 2.0)]
 [(7, [14.0, 15.0], 4.0, 4.0, 2.0)]
 [(8, [16.0, 17.0], 2.0, 2.0, 3.0)]
 [(9, [18.0, 19.0], 1.0, 2.0, 3.0)]]
|improve this answer|||||
  • I am sorry, you ar loosing me a bit. Looking at the code above, dt1 nor dt2 er defined, and final extends an object foo which isn't defined either. Maybe I am missing something but that wouldn't work as far as I can see. I would have also expected that I actually extend the distance array, rather than foo!? – Thomas Becker Jun 6 '16 at 14:09
  • Foo... or a copy of it, is what you have to append to, since cdist produces the array with reference to it. In the example I gave dist_array.shape = (10,10). final.shape is (10,) because it is a recarray and as such the columns are referenced by name. I will add an example in my post. – user681 Jun 6 '16 at 14:20
  • Thank you very much for the help Dan! Using your example I can reproduce your results. However, when I try the rec_append_fields then I am getting the following error: ValueError: The number of arrays does not match the number of names which is understandable, since I have much more start points, than I have destination points, hence my array is not square. My orig has a shape of (670911L,), dest has a shape of (39L,), resulting in a shape of dist_arr of (670911L, 39L) How can I overcome this? – Thomas Becker Jun 7 '16 at 10:44
  • adding to the link with another example – user681 Jun 7 '16 at 19:04

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.