Exclude areas (#3)

* done

Author: shrddr

lib/helpers.py

@@ -1,3 +1,4 @@
+import re
 import cv2
 import math
 import random
@@ -31,7 +32,7 @@ def outside(point, lefttop, rightbot):
         or point[1] < lefttop[1] \
         or point[0] >= rightbot[0] \
         or point[1] >= rightbot[1]
-
+    
 class MercatorPainter:
     # paints an area with dots and lines representing positive examples.
     # everything not painted over is supposed to be negative.
@@ -80,11 +81,17 @@ def add_line_wgs(self, latlng1, latlng2, width):
         p2 = self.wgs2px(latlng2)
         cv2.line(self.canvas, p1, p2, 255, width)
 
-    def add_polyline_wgs(self, latlngs, width):
-        pixels = [self.wgs2px(l) for l in latlngs]
+    def add_polyline_wgs(self, latlngs, width=1):
+        pixels = [self.wgs2px(ll) for ll in latlngs]
+        pixels = np.array(pixels)
+        # lineType means 4-connected or 8-connected
+        cv2.polylines(self.canvas, [pixels], True, 255, width, lineType=4)
+        
+    def add_fillpoly_wgs(self, latlngs):
+        pixels = [self.wgs2px(ll) for ll in latlngs]
         pixels = np.array(pixels)
         # lineType means 4-connected or 8-connected
-        cv2.polylines(self.canvas, [pixels], False, 255, width, lineType=4)
+        cv2.fillPoly(self.canvas, [pixels], 255, lineType=4)
 
     def show(self):
         # displays the canvas at native resolution
@@ -93,13 +100,14 @@ def show(self):
 
     def show_fixedwindow(self, h, w):
         # displays the canvas resized to specified dimensions
+        # useful for very small or large maps
         cv2.namedWindow('canvas-fixed', cv2.WINDOW_NORMAL)
         cv2.resizeWindow('canvas-fixed', h, w)
         cv2.imshow('canvas-fixed', self.canvas)
         cv2.waitKey(0)
 
     def build_index(self):
-        # creates lookup dict of occupied pixels
+        # creates dict of occupied pixels for fast lookup
         d = {}
         for y in range(self.height):
             ty = self.tymin + y
@@ -115,7 +123,7 @@ def build_index(self):
         self.dict_busy = d
 
     def build_index_free(self):
-        # creates lookup dict of non-occupied pixels
+        # creates dict of non-occupied pixels for even faster lookup
         d = {}
         for y in range(self.height):
             ty = self.tymin + y
@@ -205,18 +213,33 @@ def random_free(self):
         if self.dict_free[tx] == []:
             self.dict_free.pop(tx)
         return tile
+    
+def latlngs_from_wkt(string):
+    latlngs = []
+    polys = string.splitlines()
+    for poly in polys:
+        strs = re.findall(r"[-]?\d*\.\d+|\d+", poly)
+        nums = list(map(float, strs))
+        lngs = nums[::2]
+        lats = nums[1::2]
+        latlngs.append(list(zip(lats,lngs)))
+    return latlngs
 
 if __name__ == "__main__":
-    box = (27.4026,53.8306,27.7003,53.9739)
-    lamps = [(53.85, 27.6), (53.92, 27.5)]
-    roads = [[(53.85, 27.5), (53.92, 27.6)]]
+#    box = (27.4026,53.8306,27.7003,53.9739)
+#    lamps = [(53.85, 27.6), (53.92, 27.5)]
+#    roads = [[(53.85, 27.5), (53.92, 27.6)]]
+#    
+#    mp = MercatorPainter(layers.maxar, *box, z=18)
+#    mp.add_dots_wgs(lamps)  
+#    
+#    for nodes in roads:
+#        mp.add_polyline_wgs(nodes, width=2)
+#       
+#    mp.show()
+#
+#    print(mp.contains((302304, 168755)))
 
-    mp = MercatorPainter(layers.maxar, *box, z=18)
-    mp.add_dots_wgs(lamps)  
-    
-    for nodes in roads:
-        mp.add_polyline_wgs(nodes, width=2)
-       
-    mp.show()
-
-    print(mp.contains((302304, 168755)))
+    s = """POLYGON ((1.1 .2, 1 2.2, 1 -2.2))
+             POLYGON ((1 2, 1 2, 1 2))"""
+    print(latlngs_from_wkt(s))

lib/layers.py

@@ -7,6 +7,9 @@
 import numpy as np
 from pathlib import Path
 
+# every WGS coord here is (LATITUDE, LONGITUDE) - note the order!
+# LAT means up-down/north-south, LNG means left-right/west-east
+
 TILESIZE = 256
 
 def get_or_sleep(sess, url, t=0.1):
@@ -115,8 +118,8 @@ def tile_at_wgs(self, latlng, z):
         return (tx, ty)
 
     def gettile_wgs(self, latlng, z, skipedge=False, edge=16):
-        # returns tile at location (as filename string)
-        # returns None if skipedge is enabled and location is indeed close to edge 
+        # returns tile at WGS point (as filename string)
+        # returns None if skipedge is enabled and the point is indeed close to edge 
         scale = 1 << z
         wc = project2web(latlng)
         # pixel in world
@@ -126,14 +129,14 @@ def gettile_wgs(self, latlng, z, skipedge=False, edge=16):
         tx = math.floor(px / TILESIZE)
         ty = math.floor(py / TILESIZE)
         # pixel in tile
-        rx = (px - tx) % TILESIZE
-        ry = (py - ty) % TILESIZE
+        rx = px - tx * TILESIZE
+        ry = py - ty * TILESIZE
 
         if skipedge:
-            edge = (rx < edge) or (rx >= TILESIZE-edge) \
-                or (ry < edge) or (ry >= TILESIZE-edge)
-            if edge:
-                print("edge")
+            outlier = (rx < edge) or (rx >= TILESIZE-edge) \
+                   or (ry < edge) or (ry >= TILESIZE-edge)      
+            if outlier:
+                print("outlier")
                 return None
 
         fname = self.download(tx, ty, z)

make_buildings.py

@@ -14,51 +14,55 @@ def mil(fp):
     return math.floor(fp*1000000)
 
 if __name__ == "__main__":
-#    box = (27.4631,53.9359,27.4989,53.9566)
+    # bounding box to work on
     box_minsk = (27.4013,53.8157,27.7827,53.9739)
 #    box_hrodna = (23.7483,53.5909,23.9145,53.7544)   
 
+    # These areas are where active construction/destruction is 
+    # going on and the satellite imagery is outdated compared to OSM.
+    # Use http://geojson.io to draw, then export as WKT
+    with open("make_buildings_except.wkt") as reader:
+        wkt = reader.read()
+    # or just initialise with [] if you have nothing to exclude
+    outdated_polys = helpers.latlngs_from_wkt(wkt)
+          
     q = loaders.Querier()
     ways = q.query_buildings(*box_minsk)
 
     IMZ = 18
+    LIMIT = 10000
+    count = 0
     target = helpers.cleandir(f"buildings/yes")
     for wayid,nodes in ways:
         print(wayid)
-        
         # fetch every tile the building has a node at
-        for i,node in enumerate(nodes):
-            fname = layers.maxar.gettile_wgs(node, IMZ, skipedge=True, edge=16)
-            if fname is None:
+        for node in nodes:
+            fpath = layers.maxar.gettile_wgs(node, IMZ, skipedge=True, edge=24)
+            if fpath is None:
                 continue
-            dst = str(target)
-            shutil.copy(fname, dst)
-            
-        # OR focus on the building and crop
-#        image = layers.maxar.tiles_way(nodes, IMZ, pad_pct=0.25, pad_px=48)
-#        dst = str(target / f"m{wayid}.jpg")
-#        cv2.imwrite(dst, image)
+            fname = os.path.basename(fpath)
+            dst = target / fname
+            if not os.path.isfile(dst):
+                shutil.copy(fpath, dst)
+                count += 1      
+                
+        if count >= LIMIT:
+            break
 
     mp = helpers.MercatorPainter(layers.maxar, *box_minsk, IMZ)
     for _, nodes in ways:
-        mp.add_polyline_wgs(nodes, width=1)
-    # TODO: add restricted zones where imagery is outdated
-    # use cv2.polyline(is_closed=True)
-        
+        mp.add_polyline_wgs(nodes)
+    for poly in outdated_polys:
+        mp.add_fillpoly_wgs(poly)
+    
     target = helpers.cleandir(f"buildings/no")
 
     train = []
-    TRAIN = 5755
     start = time.time()
-    for i in range(TRAIN):
+    for i in range(LIMIT):
         tx, ty = mp.random_negative()
         print(tx,ty)
         fname = layers.maxar.download(tx, ty, IMZ)  
-#        wgs = layers.wgs_at_tile(tx, ty, IMZ)
-#        print(f"https://www.openstreetmap.org/edit#map=18/{wgs[0]}/{wgs[1]}")    
-#        img = cv2.imread(fname)
-#        cv2.imshow("sanity check", img)
-#        cv2.waitKey(0)
 
         dst = str(target / f"m_x{tx}y{ty}.jpg")
         if not os.path.exists(dst):

make_buildings_except.wkt

@@ -0,0 +1,7 @@
+POLYGON ((27.526330947875977 53.872418900515235, 27.523326873779297 53.869736775642515, 27.538089752197266 53.85994302920735, 27.558259963989258 53.857842263021006, 27.561392784118652 53.86892704376814, 27.552337646484375 53.87206466746561, 27.53937721252441 53.87305145207004, 27.526330947875977 53.872418900515235))
+POLYGON ((27.59085416793823 53.84343778491824, 27.588472366333008 53.842361680295205, 27.590253353118907 53.839031910806206, 27.59011387825014 53.83727196370016, 27.594963312149066 53.83596145144472, 27.598900794982928 53.83719599309501, 27.59085416793823 53.84343778491824))
+POLYGON ((27.52163171768188 53.906335460705726, 27.520623207092285 53.903908344100664, 27.52152442932129 53.90282115239138, 27.525236606597897 53.90402211834181, 27.523412704467773 53.90631017896792, 27.52163171768188 53.906335460705726))
+POLYGON ((27.55069077014923 53.90880983675384, 27.551409602165222 53.90785233697202, 27.552316188812256 53.90823786757667, 27.55150616168976 53.909008918110516, 27.55069077014923 53.90880983675384))
+POLYGON ((27.600810527801514 53.837841738845775, 27.60177612304687 53.83675282848474, 27.60791301727295 53.83719599309501, 27.609500885009766 53.83904457202817, 27.61707544326782 53.83989286518068, 27.61662483215332 53.84076646255469, 27.604007720947262 53.83910787808068, 27.600810527801514 53.837841738845775))
+POLYGON ((27.542606592178345 53.8908601761207, 27.546662092208862 53.887951598930506, 27.550975084304806 53.88825511209557, 27.54966616630554 53.88931739081768, 27.542606592178345 53.8908601761207))
+POLYGON ((27.558270692825317 53.89635430547624, 27.562283277511597 53.89519105182811, 27.56191849708557 53.89751120497534, 27.558270692825317 53.89635430547624))

make_roofshapes.py

@@ -39,17 +39,11 @@ def mil(fp):
             ways[cat] = ways[cat][:LIMIT]
         for wayid, nodes in ways[cat]:
             image = layers.maxar.tiles_way(nodes, IMZ, pad_pct=0.25, pad_px=48)
-#            print(image.shape)
             if helpers.outside(image.shape[:2], (128,128), (1024,1024)):
                 continue
-#            cv2.imshow('crop', image)
-#            cv2.waitKey(0)
-#            cv2.destroyAllWindows()
             dst = str(target / f"m{wayid}.jpg")
             cv2.imwrite(dst, image)
-#            image = layers.dg.tiles_way(nodes, IMZ)
-#            dst = str(target / f"d{wayid}.jpg")
-#            cv2.imwrite(dst, image)
+
 
     tarball = "./roofshapes.tar"
     if os.path.exists(tarball):

readme.md

@@ -1,7 +1,11 @@
-This is a collection of scripts to generate machine learning training images using OpenStreetMap data and satellite imagery providers. I tend to use [Maxar](https://github.com/osmlab/editor-layer-index/pull/655) layer because it's recent, crispy, and shot at almost vertical angle. At other locations results may vary.
+This is a collection of scripts to generate machine learning training images using OpenStreetMap data and satellite imagery providers. The scripts scan through certain objects in a specified area (defined by bounding box) and generate a .tar file with images, labeled by folder names. Upload that to your NVidia farm and train.
+
+I tend to use [Maxar](https://github.com/osmlab/editor-layer-index/pull/655) layer because it's recent, crispy, and shot at almost vertical angle. At other locations results may vary. You can add new imagery providers to `layers.py`.
 
 The training uses resnet34 architecture with fast.ai library. Every table below lists `error_rate` metric, which is just a percentage of incorrect predictions on the validation set.
 
+If you have a high error rate, fast.ai provides a great debugging tool `top_losses()`. It will give you images that confuse the net for some reason. It translates into unmapped areas or outdated imagery (or a bug in data collection script).
+
 # Streetlamps
 
 Task: Classify imagery tiles with [streetlamp](https://wiki.openstreetmap.org/wiki/Tag:highway%3Dstreet_lamp)s vs. no streetlamps. 
@@ -135,9 +139,9 @@ For every node of every building, fetch the tile it belongs to. Don't use the ti
 
 |                                  | Frozen, 4 epoch train | Unfreeze, 6 epochs | 8 more epochs |
 | -------------------------------- | --------------------- | ------------------ | ------------- |
-| size=256 no flip, max_rotate=0   | 3.0%                  | 3.0%               |               |
-| size=256 no flip, max_rotate=20  | 2.3%                  | 2.7%               |               |
+| size=256, no flip, max_rotate=0  | 3.0%                  | 3.0%               |               |
+| size=256, no flip, max_rotate=20 | 2.3%                  | 2.7%               |               |
 | size=256, no flip, max_rotate=40 | 3.3%                  | 3.3%               |               |
-|                                  |                       |                    |               |
+| size=256, no flip, max_rotate=20 | 2.5%                  |                    |               |
 |                                  |                       |                    |               |
 |                                  |                       |                    |               |