Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I have some cross section survey data, which is like this:

X Y Z pt_dt
2970969.635 359725.0088 83.4242 1-x
2970968.278 359722.2182 83.2591 1-x
2970941.771 359670.127  83.0655 1-x
2970961.369 359708.6424 83.4785 1-x
2970950.888 359688.048  85.2994 1-x
2970955.084 359696.2922 83.6479 1-x
2970953.47  359693.1165 84.9628 1-x
2970948.813 359683.977  83.1451 1-x
2970931.904 359650.7413 83.1042 1-x

When I use the X & Y to make points, I can see the points on the map like this: enter image description here

What I need to do now, is to get the points in linear order (i.e. from left to right, or right to left), so that I can create a 3D line, which can serve as the Cross sectional profile.

My current algorithm, is to calculate the distance between every pair of points(For one given cross section, which is indicated by pt_dt). The pair which has the maximum distance, are the outermost points. I then arrange the points based on distance from one of the outermost points. I have implemented this algorithm in python, and it works.

This algorithm is O(n2); My entire data consists of about 5000 such surveys, so any complexity that I can reduce will be beneficial.

So Is there a more efficient way, or algorithm to do this?

share|improve this question
3  
If your points are on an almost straight line, max/min X and Y should give you the outermost points without distance calculation of all pairs. –  AndreJ Aug 27 at 6:04
1  
@Andre A lexicographic sort would do the same but be more robust and reliable. Once the endpoints (X0,Y0) and (X1,Y1) are found, all points (X,Y) can then be sorted by their dot products with the direction vector; that is, by the values (X-X0)*(X1-X0) + (Y-Y0)*(Y1-Y0). –  whuber Aug 27 at 14:40
    
@whuber: Can you explain what 'lexicographic sort' would mean in this context? I've only ever used this sorting for strings. –  Devdatta Tengshe Aug 28 at 2:47
    
@Devdatta Sorry for being so terse: I was referring to a lexicographic sort on the x and y coordinates. (Such sorting is often a first step in many algorithms of computational geometry.) It is accomplished simply by sorting on one of the coordinates and then using the second coordinate to resolve any ties. This sort is built in to many programming environments. Finding the endpoints needs only O(n*log(n)) time and O(n) space (or O(n) time with a clever algorithm). It works well here because it automatically handles points that are lined up perfectly vertically or perfectly horizontally. –  whuber Aug 28 at 14:16

1 Answer 1

up vote 4 down vote accepted

Taking into account Andre's comment, and for the sake of learning I modified another user's answer to create a 3D Line from your sample data. I hope you or someone else might find it useful:

enter image description here

X,Y,Z,pt_dt
2970969.635,359725.0088,83.4242,1-x
2970968.278,359722.2182,83.2591,1-x
2970941.771,359670.127,83.0655,1-x
2970961.369,359708.6424,83.4785,1-x
2970950.888,359688.048,85.2994,1-x
2970955.084,359696.2922,83.6479,1-x
2970953.47,359693.1165,84.9628,1-x
2970948.813,359683.977,83.1451,1-x
2970931.904,359650.7413,83.1042,1-x

import csv
import numpy as np
try:
    from osgeo import ogr,osr
except ImportError:
    import ogr,osr

SHP_FILENAME = "output.shp"
CSV_FILENAME = "exp_data.csv"

r = csv.reader(open(CSV_FILENAME, 'r'), delimiter=',', quotechar=None)
header = dict(((str, i) for i, str in enumerate(r.next())))
#print header

# load data rows into memory
rows = [row for row in r]


# sort by x if larger than y. otherwise sort by y
np_rows = np.array(rows)
np_rows = np.delete(np_rows,header['pt_dt'],axis=1)
np_rows = np_rows.astype(np.float)

rows.sort(key=lambda x: x[0] if np_rows.ptp(axis=0)[0] > np_rows.ptp(axis=0)[1] else x[1])


# Create new shapefile
ogr.UseExceptions()
sr = osr.SpatialReference()
sr.ImportFromEPSG(900913)
ds = ogr.GetDriverByName('ESRI Shapefile').CreateDataSource(SHP_FILENAME)

layer = ds.CreateLayer("data", sr, ogr.wkbLineString25D)

# Create a new line geometry
line = ogr.Geometry(type=ogr.wkbLineString25D)
line.SetCoordinateDimension(3)



# Add GPS points to line
lon_idx, lat_idx ,Z_idx = header['X'], header['Y'], header['Z']
for row in rows:
    line.AddPoint(float(row[lon_idx]),float(row[lat_idx]),float(row[Z_idx]))

# Add line as a new feature to the shapefile
feature = ogr.Feature(feature_def=layer.GetLayerDefn())
feature.SetGeometryDirectly(line)

layer.CreateFeature(feature)

# Cleanup
feature.Destroy()
ds.Destroy()
share|improve this answer
    
If the cross section line is North-South with some tolerance, your sorting will fail. You better look at maxX-minX and maxY-minY, and take the axis of the larger extent for sorting. –  AndreJ Aug 27 at 9:28
    
Good point, I've updated the script –  nickves Aug 27 at 10:47

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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