0

I have a point layer with site_code attribute data. I have another layer of points in different locations but with the same site_code data.

How can I find the linear distance between points that are paired by their site_code attribute data?

  • 3
    What software/version are you using? Also, are you looking for GUI or script answer? – phloem Jan 28 '15 at 17:53
3

ArcGIS GUI solution:

  1. Merge all points into one file
  2. Run Points to Line tool, using site_code as Line Field. If you save it to a geodatabase, SHAPE_LENGTH will be the distance of the line. Otherwise, calculate the length of the line yourself.

ArcGIS (10.1+) field calculator solution:

This matches points with the same 'site_code' between two files.

Parser:

Python

Expression:

findDistance( !site_code!, !Shape! )

Code Block:

def findDistance(sc, shape):
    fc = r'C:\junk\points_copy.shp' # change to your 'match file' path
    with arcpy.da.SearchCursor(fc,['site_code','SHAPE@']) as cursor:
        for row in cursor:
            if sc == row[0]:
                return row[1].distanceTo(shape)
2

With Fiona and shapely (in pure Python, no need of ArcPy or PyQGIS)

from shapely import shape
import fiona
# reading the 2 shapefiles
feat1 = list(fiona.open('point1.shp')) # list of dictionaries
feat2 = list(fiona.open('point2.shp'))
# find the linear distance between points that are paired by their site_code attribute data
pairs = zip(feat1,feat2)
print [(x['properties']['site-code'], shape(x['geometry']).distance(shape(y['geometry']))) for x, y in pair if x['properties']['site-code']== y['properties']['site-code']]
[(5, 1806.8754790315722), (4, 500.70397719211178), (3, 543.37340530091308), (2, 584.88140088639977), (1, 291.35637024434374)]

Explication:

c = fiona.open('point1.shp')
# schema of the shapefile
c.schema
{'geometry': 'Point', 'properties': OrderedDict([(u'id', 'int:10'), (u'site-code', 'int:5')])}
# first feature of the shapefile
first = c.next()
print first
{'geometry': {'type': 'Point', 'coordinates': (151298.25462051606, 132341.01651859211)}, 'type': 'Feature', 'id': '0', 'properties': OrderedDict([(u'id', None), (u'site-code', 5)])}
# the result is a dictionary, so:`
print first['properties']['site-code']
5
print shape(first['geometry'])
POINT (151298.2546205161 132341.0165185921) # a shapely geometry

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