organizing, unit testing, rewriting, likely also breaking working stu…

…ff. #5

ghini_tree_position.py

@@ -207,92 +207,46 @@ 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()
             # force editable if not already editable
             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)

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)