research

Author: shrddr

helpers.py

@@ -19,13 +19,14 @@ class MercatorPainter:
     # everything not painted over is supposed to be negative.
     # uses dict for fast lookup (builds itself on first query).
     # also has a function to find a random negative (unpainted) pixel.
-    def __init__(self, W, S, E, N):   
-        txmin, tymin = layers.tile_at_wgs((N, W))
-        txmax, tymax = layers.tile_at_wgs((S, E))
+    def __init__(self, layer, W, S, E, N):   
+        txmin, tymin = layer.tile_at_wgs((N, W))
+        txmax, tymax = layer.tile_at_wgs((S, E))
         area = (txmax-txmin, tymax-tymin)
         print(f"paint area: {txmin}..{txmax}, {tymin}..{tymax}")
         print(f"dimensions: {area} -> {area[0]*area[1]} tiles total")
 
+        self.layer = layer
         self.txmin = txmin
         self.tymin = tymin
         self.width = txmax-txmin+1
@@ -35,7 +36,7 @@ def __init__(self, W, S, E, N):
         self.dict = None
 
     def wgs2px(self, latlng):
-        tx, ty = layers.tile_at_wgs(latlng)
+        tx, ty = self.layer.tile_at_wgs(latlng)
         x = tx - self.txmin
         y = ty - self.tymin
         return (x,y)
@@ -46,18 +47,18 @@ def add_dot_tile(self, tile, color=255):
         y = ty - self.tymin
         self.canvas[y][x] = color
 
-    def add_dots(self, latlngs, color=255):
+    def add_dots_wgs(self, latlngs, color=255):
         for latlng in latlngs:
             x, y = self.wgs2px(latlng)
             self.canvas[y][x] = color
 
-    def add_line(self, latlng1, latlng2, width):
+    def add_line_wgs(self, latlng1, latlng2, width):
         # the lines are actually curves in mercator but we don't care
         p1 = self.wgs2px(latlng1)
         p2 = self.wgs2px(latlng2)
         cv2.line(self.canvas, p1, p2, 255, width)
 
-    def add_polyline(self, latlngs, width):
+    def add_polyline_wgs(self, latlngs, width):
         pixels = [self.wgs2px(l) for l in latlngs]
         pixels = np.array(pixels)
         cv2.polylines(self.canvas, [pixels], False, 255, width)
@@ -108,6 +109,7 @@ def find_random(self):
             if self.contains(tile):
                 continue
             else:
+                self.add_dot_tile(tile)
                 return tile
 
 if __name__ == "__main__":

layers.py

@@ -16,14 +16,6 @@ def project2web(latlng):
     y = TILESIZE * (0.5 - math.log((1 + siny) / (1 - siny)) / (4 * math.pi))
     return (x, y)
 
-def tile_at_wgs(latlng, z=19):
-    # returns index of tile which contains a location
-    scale = 1 << z
-    wc = project2web(latlng)
-    tx = math.floor(wc[0] * scale / TILESIZE)
-    ty = math.floor(wc[1] * scale / TILESIZE)
-    return (tx,ty)
-
 def wgs_at_tile(tx, ty, z=19):
     # converts tile index to EPSG:3857 (0..1) then to EPSG:4326 (degrees)
     scale = 1 << z
@@ -41,7 +33,7 @@ def __init__(self, name):
         self.flipy = False
         self.offsetx = 0
         self.offsety = 0
-        self.tiledir = Path("tiles") / name
+        self.tiledir = Path("../tiles") / name
 
     def tilefile(self, x, y, z):
         return self.tiledir / f"x{x}y{y}z{z}.jpg"
@@ -72,10 +64,42 @@ def download(self, x, y, z=19):
                 raise IOError(f"{r.status_code} at {url}'")
         return str(fname)
 
-    def gettile_wgs(self, latlng, z=19):
+    def tile_at_wgs(self, latlng, z=19):
+        # returns index of tile which contains a location
+        scale = 1 << z
+        wc = project2web(latlng)
+        # pixel in world
+        px = wc[0] * scale + self.offsetx
+        py = wc[1] * scale + self.offsety
+        # tile in world
+        tx = math.floor(px / TILESIZE)
+        ty = math.floor(py / TILESIZE)
+        return (tx, ty)
+    
+    def gettile_wgs(self, latlng, z=19, skipedge=False):
         # returns tile at location (as filename)
-        x,y = tile_at_wgs(latlng, z)
-        fname = self.download(x, y, z)
+        # returns None if skipedge is enabled and location is indeed close to edge 
+        scale = 1 << z
+        wc = project2web(latlng)
+        # pixel in world
+        px = wc[0] * scale + self.offsetx
+        py = wc[1] * scale + self.offsety
+        # tile in world
+        tx = math.floor(px / TILESIZE)
+        ty = math.floor(py / TILESIZE)
+        # pixel in tile
+        rx = (px - tx) % TILESIZE
+        ry = (py - ty) % TILESIZE
+
+        if skipedge:
+            EDGE_THRESH = 10 # px
+            edge = (rx < EDGE_THRESH) or (rx >= TILESIZE-EDGE_THRESH) \
+                or (ry < EDGE_THRESH) or (ry >= TILESIZE-EDGE_THRESH)
+            if edge:
+                print("edge")
+                return None
+
+        fname = self.download(tx, ty, z)
         return fname
 
     def tiles_near_wgs(self, latlng, scale, h, w):
@@ -148,4 +172,4 @@ def getcrop_wgs(self, latlng, h, w, z=19):
 dg = Imagery("dg")
 dg.url = "https://c.tiles.mapbox.com/v4/digitalglobe.316c9a2e/{z}/{x}/{y}.png?access_token=pk.eyJ1IjoiZGlnaXRhbGdsb2JlIiwiYSI6ImNqZGFrZ2c2dzFlMWgyd2x0ZHdmMDB6NzYifQ.9Pl3XOO82ArX94fHV289Pg"
 
-print(maxar.xy_fromfile(Path(r"tiles\maxar\x302117y168688z19.jpg")))
+print(maxar.xy_fromfile(Path(r"..\tiles\maxar\x302117y168688z19.jpg")))

loaders.py

@@ -1,26 +1,21 @@
 import overpass
 import json
 import os.path
-import time
-import datetime
 from math import floor
 
-import helpers
-import layers
-
 def mil(fp):
     return floor(fp*1000000)
 
 def query_nodes(W, S, E, N):
     # queries overpass or fetches cached result if available
     # returns list of (lat, lng) tuples
-    fname = f"./in/bbox{mil(W)}_{mil(S)}_{mil(E)}_{mil(N)}.json"
+    fname = f"./overpass/bbox{mil(W)}_{mil(S)}_{mil(E)}_{mil(N)}.json"
     if os.path.isfile(fname):
         with open(fname) as json_file:
             return json.load(json_file)
 
     api = overpass.API()
-    query = f"node[\"highway\"=\"street_lamp\"]({S}, {W}, {N}, {E})"
+    query = f"""node["highway"="street_lamp"]({S}, {W}, {N}, {E})"""
     response = api.get(query, responseformat='json', verbosity='skel')
 
     nodes = [(e['lat'], e['lon']) for e in response['elements']]
@@ -32,7 +27,7 @@ def query_nodes(W, S, E, N):
 
 
 def query_ways(W, S, E, N):
-    fname = f"./in/ways_bbox{mil(W)}_{mil(S)}_{mil(E)}_{mil(N)}.json"
+    fname = f"./overpass/ways_bbox{mil(W)}_{mil(S)}_{mil(E)}_{mil(N)}.json"
     if os.path.isfile(fname):
         with open(fname) as json_file:
             return json.load(json_file)
@@ -74,6 +69,5 @@ def query_poly(bounds):
     # https://wiki.openstreetmap.org/wiki/Overpass_API/Overpass_QL#Polygon_evaluator
     pass
 
-
 if __name__ == "__main__":   
     pass

make.py

@@ -10,7 +10,10 @@
 import tarfile
 
 # use one original tile or expand it for later cropping
-ORIGINAL = False
+MAKE_POSITIVE_ORIGINAL = 0
+MAKE_POSITIVE_EXPANDED = 0
+MAKE_NEGATIVE_ORIGINAL = 0
+MAKE_NEGATIVE_EXPANDED = 1
 
 # this is 256 for all current imagery providers
 TILESIZE = 256
@@ -19,7 +22,7 @@
 # known satellite imagery has up to 40px offset.
 # maximum acceptable padding is (128-40)=88
 # which translates to image size 256+88*2=432
-PADDING = 50
+PADDING = 88
 
 if __name__ == "__main__": 
 #    box = (27.4583,53.9621,27.5956,53.9739) # north belt
@@ -31,19 +34,20 @@
     lamps = loaders.query_nodes(*box)
     print("lamps:", len(lamps))
 
-    if ORIGINAL: 
-        # only use one tile where lamp exists. should fail
-        # miserably when the lamp is on the edge of tile
+    if MAKE_POSITIVE_ORIGINAL: 
+        # only use one tile where lamp exists
         target = pathlib.Path('lamps-orig/lamp')
         target.mkdir(parents=True, exist_ok=True)
 
         for lamp in lamps:
-            fname = layers.maxar.gettile_wgs(lamp)
-            dst = target / ("m_" + os.path.basename(fname))
-            shutil.copy(fname, dst)
+            fname = layers.maxar.gettile_wgs(lamp, skipedge=True)
+            if fname is not None:
+                dst = target / ("m_" + os.path.basename(fname))
+                shutil.copy(fname, dst)
 
-    else:
-        # use a bigger picture for later random cropping (augmentation)
+    if MAKE_POSITIVE_EXPANDED:
+        # use a bigger picture, centered at the object, 
+        # for later random cropping (augmentation)
         target = pathlib.Path('lamps-center/lamp')
         target.mkdir(parents=True, exist_ok=True)
 
@@ -61,11 +65,11 @@
 
     BATCHSIZE = len(lamps)
 
-    mp = helpers.MercatorPainter(*box)
-    mp.add_dots(lamps)
+    mp = helpers.MercatorPainter(layers.maxar, *box)
+    mp.add_dots_wgs(lamps)
     roads = loaders.query_ways(*box)
     for nodes in roads.values():
-        mp.add_polyline(nodes, width=2)
+        mp.add_polyline_wgs(nodes, width=2)
 
     source = layers.maxar.tiledir
     localtiles = [layers.maxar.xy_fromfile(path) for path in source.glob("*.jpg")]
@@ -82,38 +86,37 @@
         while BATCHSIZE > len(batch):
             batch.append(mp.find_random())
 
-    if ORIGINAL: # only use one tile
-        # TODO: discard if lamp is close to the edge
-        # TODO: handle imagery offset
+    if MAKE_NEGATIVE_ORIGINAL: 
+        # only use one tile
         target = pathlib.Path('lamps-orig/nolamp')
         target.mkdir(parents=True, exist_ok=True) 
         for (tx,ty) in batch:
             fname = layers.maxar.download(tx,ty)
             dst = target / ("m_" + os.path.basename(fname))
             shutil.copy(fname, dst)
 
-    else:
+    if MAKE_NEGATIVE_EXPANDED:
         # expand the tile for later cropping
         target = pathlib.Path('lamps-center/nolamp')
         target.mkdir(parents=True, exist_ok=True) 
         for (tx,ty) in batch:
             wgs = layers.wgs_at_tile(tx, ty)
-            h, w = (356, 356)
+            h = w = PADDING + TILESIZE + PADDING
             crop = layers.maxar.getcrop_wgs(wgs, h, w)
             lat = helpers.mil(wgs[0])
             lng = helpers.mil(wgs[1])
             dst = str(target / f"m_lat{lat}lng{lng}z19.jpg")
             cv2.imwrite(dst, crop)
 
-    if ORIGINAL:
+    if MAKE_POSITIVE_ORIGINAL and MAKE_NEGATIVE_ORIGINAL:
         tarball = "./lamps-orig.tar"
         if os.path.exists(tarball):
             os.remove(tarball)
         tar = tarfile.open(tarball, "w")
         tar.add("./lamps-orig")
         tar.close()
 
-    else:
+    if MAKE_POSITIVE_EXPANDED and MAKE_NEGATIVE_EXPANDED:
         tarball = "./lamps-center.tar"
         if os.path.exists(tarball):
             os.remove(tarball)

in/.gitignore renamed to overpass/.gitignore

@@ -1,32 +1,68 @@
+## What's up?
+
+These scripts prepare the data to teach a classification network which tells apart satellite imagery tiles with [streetlamp](https://wiki.openstreetmap.org/wiki/Tag:highway%3Dstreet_lamp)s or with no streetlamps. The tests only use the latest and the clearest imagery layer for my city called Maxar.
+
 ## Original tiles
 
-These use skew and rotate to augment (which is hardly acceptable for satellite imagery). Also, must be problematic when the lamp is on tile edge. Sometimes imagery offset makes the object appear on a wrong tile, this is probably a major roadblock.
+### First attempt
+
+Satellite imagery providers serve data in 256x256 tiles. The first approach is just fetch a tile which contains a lamp and use that as a positive example. Every tile that does not contain a high-level road (highway=tertiary and up) is supposed to be negative.
+
+The problem is some lamps are at the tile edge and possibly cross the boundary. Sometimes [imagery offset](https://wiki.openstreetmap.org/wiki/Using_Imagery#Frequent_mistakes) makes the object appear on a different tile than it should, which produces false positive example.
 
 ```
 tfms = get_transforms(do_flip=False)
 data = ImageDataBunch.from_folder(path, train=".", valid_pct=0.1, ds_tfms=tfms, size=256)
 ```
 
-Best I got was 3.9% error
+Best I got was 3.1-3.9% validation error and it just doesn't train any further.
+
+### Discard "edge cases"
+
+Here I drop all positive examples where the base of street lamp is less than 10px away from tile edge. Dealing with offset is tricky, as it depends on both imagery properties and OSM mappers in the area. I did eyeball an average offset for my area, but for other cities it can be anything. You can use `video.py` to do that, it lets you look through lots of imagery quickly.
+
+Negative examples are still one random non-road tile.
+
+This dataset converges to 3% error.
 
 ## Expanded tiles
 
-This is where I turn off skew, rotate and zoom, because crop should be sufficient. Turns out crop doesn't work like that and must be turned on explicitly - see below. This however shows the performance with virtually no augmentation.
+### Damn it works.. oh wait
+
+The best thing about satellite imagery is that it's huge and can be scrolled in every direction almost infinitely. If you need more information about a location you can always look at adjacent tiles. This method just fetches a larger square of 356 pixels around every known streetlamp. These will get randomly cropped later in training process. Negative examples are expanded to 356px too, just for consistency.
+
+Turns out, these runs had no crop at all - see below. This however shows the performance with virtually no augmentation.
 
 ```
-tfms = get_transforms(do_flip=False, max_warp=0, max_zoom=0, max_rotate=0
+# just resizes 356->256px with no crop!
+tfms = get_transforms(do_flip=False, max_warp=0, max_zoom=0, max_rotate=0)
 data = ImageDataBunch.from_folder(path, train=".", valid_pct=0.1, ds_tfms=tfms, size=256)
 ```
 
-432px: 1.3% error
+The validation error comes down to:
 
-356px: frozen 1ep 3.8% > unfrozen 2ep 1.3-1.7% > 2ep 0.7-1.0%
+|       | Deeper layers frozen, 1 epoch train | Unfreeze, 2 epochs | 2 more epochs |
+| ----- | ----------------------------------- | ------------------ | ------------- |
+| 432px |                                     |                    | 1.3%          |
+| 356px | 3.8%                                | 1.3-1.7%           | 0.7-1.0%      |
+
+### The real thing
+
+Train set cropping seems important for real-life applications, because input at inference time will include streetlamps at any part of a tile, not just center of it. That's why we need to crop randomly at train time. Random cropping turned out tricky in fast.ai, here's my hack that looks like it's doing the right thing:
 
 ```
-tfms = [[crop(size=256, row_pct=(0,1), col_pct=(0,1))],[]]
+# crops random 256x256 piece out of larger 
+# input image both at train and validation
+tfms = [[crop(size=256, row_pct=(0,1), col_pct=(0,1))],
+        [crop(size=256, row_pct=(0,1), col_pct=(0,1))]]
 data = ImageDataBunch.from_folder(path, train=".", valid_pct=0.1, ds_tfms=tfms, size=None)
 ```
 
-356px: frozen 1ep 4-6% > unfrozen 2ep 0.1-0.4% > 2ep 0.2%
+Results are indeed much better:
+
+|       | Deeper layers frozen, 1 epoch train | Unfreeze, 2 epochs | 2 more epochs |
+| ----- | ----------------------------------- | ------------------ | ------------- |
+| 356px | 4-6%                                | 0.1-0.4% / 3.4%    | 0.2%          |
+| 432px | 4.8-5.5%                            | 3.8-4.7%           | 5.5%          |
 
-As my set only contains 10k images, this means just 2 incorrect validation results O_o
+As my set only contains 10k images, and validation is 10% of that, 0.1% error means just one incorrect validation result. Oh my.