diff --git a/ghini_tree_position.py b/ghini_tree_position.py index 253f906..4be6eea 100644 --- a/ghini_tree_position.py +++ b/ghini_tree_position.py @@ -207,72 +207,32 @@ def run(self): if result: # get reference points from the layer layer = layers[self.dlg.comboBox.currentIndex()] + + # use position of first feature in layer to select local utm + f1 = layer.getFeatures().next() + + local_utm = QgsCoordinateReferenceSystem() + local_utm.createFromProj4(utm_zone_proj4(f1.geometry().asPoint())) + transf = QgsCoordinateTransform(layer.crs(), local_utm) + back_transf = QgsCoordinateTransform(local_utm, layer.crs()) + + # populate 'points' dict with projected coordinates points = {} for feature in layer.getFeatures(): - # TODO target coordinates system should come from layer - transf = QgsCoordinateTransform( - QgsCoordinateReferenceSystem(4326), - QgsCoordinateReferenceSystem(3117)) + # we work with local utm projection easting_northing = transf.transform( feature.geometry().asPoint()) point_id = feature['id'] points[point_id] = {'id': point_id, 'coordinates': easting_northing} - # the name of the file comes from the dialog box - distances = {} + # get distances from csv file, and compute connectivity to + # referenced points with open(self.dlg.lineEdit.text()) as f: - for l in f.readlines(): - l = l.strip() - try: - from_id, to_id, distance = l.split(',')[:3] - distance = float(distance) - except Exception, e: - print '»', l, '«', type(e), e - continue - distances.setdefault(from_id, {}) - distances.setdefault(to_id, {}) - distances[from_id][to_id] = distance - distances[to_id][from_id] = distance - points.setdefault(to_id, {'id': to_id, - "type": "Point"}) - points.setdefault(from_id, {'id': from_id, - "type": "Point"}) - - # inform each point on how many links lead to referenced point - for n, point in points.items(): - point['prio'] = len( - filter(lambda x: 'coordinates' in points[x], - distances.get(n, {}).keys())) - point['computed'] = False - - # remember last attempted point, to avoid deadlocks - last_attempted_point = None - # construct priority queue of points for which we still have no - # coordinates - heap = Heap([p for p in points.values() if 'coordinates' not in p]) - - while heap: - point = heap.pop() - # compute coordinates of point - try: - point['coordinates'] = list( - find_point_coordinates(points, distances, point['id'])) - except ValueError: - point['prio'] = 2 - if last_attempted_point != point: - heap.push(point) - last_attempted_point = point - continue - point['computed'] = True - - # inform points connected to point that they have one more - # referenced neighbour - for neighbour_id, destinations in distances[ - point['id']].items(): - neighbour = points[neighbour_id] - if 'heappos' in neighbour: - heap.reprioritize(neighbour) + distances = get_distances_from_csv(f, points) + + # compute missing coordinates + extrapolate_coordinates(points, distances) # remember editable status wasEditable = layer.isEditable() @@ -280,19 +240,13 @@ def run(self): if not wasEditable: layer.startEditing() - # TODO source coordinate reference system should be from active - # layer - transf = QgsCoordinateTransform( - QgsCoordinateReferenceSystem(3117), - QgsCoordinateReferenceSystem(4326)) - fields = layer.fields() featureList = [] # now add the computed points to the layer for p in [p for p in points.values() if p['computed']]: x, y = p['coordinates'] - layerPoint = transf.transform(QgsPoint(x, y)) + layerPoint = back_transf.transform(QgsPoint(x, y)) feature = QgsFeature(fields) feature.setGeometry(QgsGeometry.fromPoint(layerPoint)) feature['id'] = p['id'] @@ -320,6 +274,68 @@ def run(self): layer.commitChanges() +def extrapolate_coordinates(points, distances): + """compute missing coordinates respecting distances and given points + + navigate distances graph, keep selecting most connected point, to + compute its coordinates given enough distances from enough referenced + points. + + """ + # remember last attempted point, to avoid deadlocks + last_attempted_point = None + # construct priority queue of points for which we still have no + # coordinates + heap = Heap([p for p in points.values() if 'coordinates' not in p]) + + while heap: + point = heap.pop() + # compute coordinates of point + try: + point['coordinates'] = list( + find_point_coordinates(points, distances, point['id'])) + except ValueError: + point['prio'] = 2 + if last_attempted_point != point: + heap.push(point) + last_attempted_point = point + continue + point['computed'] = True + + # inform points connected to point that they have one more + # referenced neighbour + for neighbour_id, destinations in distances[point['id']].items(): + neighbour = points[neighbour_id] + if 'heappos' in neighbour: + heap.reprioritize(neighbour) + + +def get_distances_from_csv(stream, points): + distances = {} + for l in stream.readlines(): + l = l.strip() + try: + from_id, to_id, distance = l.split(',')[:3] + distance = float(distance) + except Exception, e: + print '»', l, '«', type(e), e + continue + distances.setdefault(from_id, {}) + distances.setdefault(to_id, {}) + distances[from_id][to_id] = distance + distances[to_id][from_id] = distance + points.setdefault(to_id, {'id': to_id, + "type": "Point"}) + points.setdefault(from_id, {'id': from_id, + "type": "Point"}) + # inform each point on how many links lead to referenced point + for n, point in points.items(): + point['prio'] = len(filter(lambda x: 'coordinates' in points[x], + distances.get(n, {}).keys())) + point['computed'] = False + return distances + + class Heap: def __init__(self, elems): """creates a binary heap from list of dictionaries @@ -608,3 +624,18 @@ def target(x): import scipy.optimize optres = scipy.optimize.minimize(target, (0, 0, 0)) return optres.x + + +def utm_zone_proj4(pt): + import math + lon, lat = pt + wkt = '+proj=utm +ellps=WGS84 +datum=WGS84 +units=m +no_defs' + zone_number = int(math.floor((lon + 180) % 360 / 6) + 1) + wkt += " +zone=%s" % zone_number + if lat < 0: + wkt += ' +south' + return wkt + + +def solve_puzzle(points, distances, gps): + place_initial_three_points(points, distances, gps) diff --git a/test/test_solve_puzzle.py b/test/test_solve_puzzle.py new file mode 100644 index 0000000..9d55f9a --- /dev/null +++ b/test/test_solve_puzzle.py @@ -0,0 +1,159 @@ +# coding=utf-8 + +import unittest +from numpy.testing import assert_almost_equal +from ghini_tree_position import solve_puzzle, get_distances_from_csv +from ghini_tree_position import utm_zone_proj4 +from ghini_tree_position import extrapolate_coordinates + + +class TestUTMChoice(unittest.TestCase): + + def test_utm_zone_proj4_europe(self): + # Naples + ex = '+proj=utm +ellps=WGS84 +datum=WGS84 +units=m +no_defs +zone=33' + self.assertEquals(utm_zone_proj4((14, 41)), ex) + # Netherlands + ex = '+proj=utm +ellps=WGS84 +datum=WGS84 +units=m +no_defs +zone=31' + self.assertEquals(utm_zone_proj4((5, 52)), ex) + # Małopolska + ex = '+proj=utm +ellps=WGS84 +datum=WGS84 +units=m +no_defs +zone=34' + self.assertEquals(utm_zone_proj4((20, 50)), ex) + + def test_utm_zone_proj4_morocco(self): + # مراكش + ex = '+proj=utm +ellps=WGS84 +datum=WGS84 +units=m +no_defs +zone=29' + self.assertEquals(utm_zone_proj4((-8, 31)), ex) + # فاس + ex = '+proj=utm +ellps=WGS84 +datum=WGS84 +units=m +no_defs +zone=30' + self.assertEquals(utm_zone_proj4((-5, 34)), ex) + + def test_utm_zone_proj4_colombia_venezuela_ecuador_panama(self): + # Quito + ex = '+no_defs +zone=17 +south' + self.assertTrue(utm_zone_proj4((-78.5, -0.2)).endswith(ex)) + # Cartagena + ex = '+no_defs +zone=18' + self.assertTrue(utm_zone_proj4((-75, 10)).endswith(ex)) + # David + ex = '+no_defs +zone=17' + self.assertTrue(utm_zone_proj4((-82.4, 8.5)).endswith(ex)) + # Iquitos + ex = '+no_defs +zone=18 +south' + self.assertTrue(utm_zone_proj4((-73, -3.7)).endswith(ex)) + + +class TestGetDistancesFromCsv(unittest.TestCase): + + def test_get_distances_from_csv(self): + points = {} + from StringIO import StringIO + s = StringIO('14,24,3\n14,25,4\n24,25,5\n') + d = get_distances_from_csv(s, points) + self.assertTrue('14' in points) + self.assertTrue('24' in points) + self.assertTrue('25' in points) + self.assertEquals(d['14']['24'], 3) + self.assertEquals(d['24']['14'], 3) + self.assertEquals(d['14']['25'], 4) + self.assertEquals(d['25']['14'], 4) + self.assertEquals(d['24']['25'], 5) + self.assertEquals(d['25']['24'], 5) + + def test_get_distances_from_csv_adding(self): + points = {'15': {'coordinates': (0, 0)}} + from StringIO import StringIO + s = StringIO('14,15,3\n') + d = get_distances_from_csv(s, points) + self.assertTrue('14' in points) + self.assertTrue('15' in points) + self.assertFalse('coordinates' in points['14']) + self.assertTrue('coordinates' in points['15']) + self.assertEquals(d['14']['15'], 3) + self.assertEquals(d['15']['14'], 3) + + def test_get_distances_connectivity_zero(self): + points = {} + from StringIO import StringIO + s = StringIO('14,24,3\n14,25,4\n24,25,5\n') + d = get_distances_from_csv(s, points) + self.assertEquals(points['14']['prio'], 0) + self.assertEquals(points['24']['prio'], 0) + self.assertEquals(points['25']['prio'], 0) + + def test_get_distances_connectivity_three(self): + points = {'15': {'coordinates': (0, 0)}, + '16': {'coordinates': (0, 0)}, + '17': {'coordinates': (0, 0)}} + from StringIO import StringIO + s = StringIO('14,15,3\n14,16,3\n14,17,3\n15,16,3\n') + d = get_distances_from_csv(s, points) + self.assertEquals(points['14']['prio'], 3) + self.assertEquals(points['15']['prio'], 1) + self.assertEquals(points['17']['prio'], 0) + + def test_get_distances_connectivity_one(self): + points = {'15': {'coordinates': (0, 0)}} + from StringIO import StringIO + s = StringIO('14,15,3\n') + d = get_distances_from_csv(s, points) + self.assertEquals(points['14']['prio'], 1) + self.assertEquals(points['15']['prio'], 0) + + +class TestExtrapolateCoordinates(unittest.TestCase): + + def test_extrapolate_coordinates_one(self): + points = {'0': {'coordinates': (0, 0)}, + 'A': {'coordinates': (4, 0)}, + 'B': {'coordinates': (0, 3)}} + from StringIO import StringIO + s = StringIO('x,0,5\nx,A,3\nx,B,4\n') + d = get_distances_from_csv(s, points) + extrapolate_coordinates(points, d) + self.assertEquals(tuple(points['x']['coordinates']), (4.0, 3.0)) + self.assertEquals(points['x']['computed'], True) + + def test_extrapolate_coordinates_two(self): + points = {'0': {'coordinates': (4, 0)}, + 'A': {'coordinates': (0, 0)}, + 'B': {'coordinates': (0, 3)}, + 'C': {'coordinates': (0, 6)}, + 'D': {'coordinates': (0, 9)}} + from StringIO import StringIO + s = StringIO('x,0,3\nx,A,5\nx,B,4\n' + 'y,x,3\ny,B,5\ny,C,4\n' + 'z,y,3\nz,C,5\nz,D,4\n' + 't,x,2.5\nt,B,2.5\nt,y,2.5\n') + d = get_distances_from_csv(s, points) + extrapolate_coordinates(points, d) + self.assertEquals(tuple(points['x']['coordinates']), (4.0, 3.0)) + self.assertEquals(tuple(points['y']['coordinates']), (4.0, 6.0)) + self.assertEquals(tuple(points['z']['coordinates']), (4.0, 9.0)) + assert_almost_equal(tuple(points['t']['coordinates']), (2.0, 4.5)) + + +class TestSolvePuzzle(unittest.TestCase): + + def test_solve_puzzle(self): + points = { + 'p1': {'id': 'p1'}, + 'p2': {'id': 'p2'}, + 'p3': {'id': 'p3'}, + 'p4': {'id': 'p4'}} + proj_gps = {'p1': {'coordinates': (0.0, 0.0)}, + 'p2': {'coordinates': (4.0, 0.0)}, + 'p3': {'coordinates': (4.0, 3.0)}, + 'p4': {'coordinates': (0.0, 3.0)}, } + distances_list = [('p1', 'p2', 4.0), + ('p1', 'p3', 5.0), + ('p1', 'p4', 3.0), + ('p2', 'p3', 3.0), + ('p2', 'p4', 5.0), + ('p1', 'p4', 4.0)] + distances = {} + for p, q, d in distances_list: + distances.setdefault(p, {}) + distances.setdefault(q, {}) + distances[p][q] = distances[q][p] = d + solution = solve_puzzle(points, distances, proj_gps)