From e0468cd79c239f62b8358c2b32ed145ead6e1269 Mon Sep 17 00:00:00 2001
From: Ian Rees <ian@ianrees.net>
Date: Tue, 27 Nov 2018 16:13:49 -0800
Subject: [PATCH 01/13] GeoTIFF downloads (WIP)

---
 planetutils/download.py                  |  8 +-
 planetutils/elevation_tile_download.py   |  4 +-
 planetutils/elevation_tile_downloader.py | 99 ++++++++++++++++--------
 3 files changed, 73 insertions(+), 38 deletions(-)

diff --git a/planetutils/download.py b/planetutils/download.py
index c586aca..2b1e033 100644
--- a/planetutils/download.py
+++ b/planetutils/download.py
@@ -1,11 +1,15 @@
 from __future__ import absolute_import, unicode_literals
 import os
 import subprocess
+import requests
 from . import log
 
 def download(url, outpath):
-    pass
-
+    r = requests.get(url, stream=True)
+    with open(outpath, 'wb') as fd:
+        for chunk in r.iter_content(chunk_size=128):
+            fd.write(chunk)
+    
 def download_gzip(url, outpath):
     with open(outpath, 'wb') as f:
         ps1 = subprocess.Popen(['curl', '-L', '--fail', '-s', url], stdout=subprocess.PIPE)
diff --git a/planetutils/elevation_tile_download.py b/planetutils/elevation_tile_download.py
index 0728e80..e779ca3 100644
--- a/planetutils/elevation_tile_download.py
+++ b/planetutils/elevation_tile_download.py
@@ -4,7 +4,7 @@
 
 from . import log
 from .bbox import load_bboxes_csv, bbox_string
-from .elevation_tile_downloader import ElevationTileDownloader
+from .elevation_tile_downloader import ElevationGeotiffDownloader
 
 def main():
     parser = argparse.ArgumentParser()
@@ -17,7 +17,7 @@ def main():
     if args.verbose:
         log.set_verbose()
 
-    p = ElevationTileDownloader(args.outpath)
+    p = ElevationGeotiffDownloader(args.outpath)
     if args.csv:
         p.download_bboxes(load_bboxes_csv(args.csv))
     elif args.bbox:
diff --git a/planetutils/elevation_tile_downloader.py b/planetutils/elevation_tile_downloader.py
index fa243f8..60c0aa3 100644
--- a/planetutils/elevation_tile_downloader.py
+++ b/planetutils/elevation_tile_downloader.py
@@ -14,9 +14,7 @@ def makedirs(path):
     except OSError as e:
         pass
 
-class ElevationTileDownloader(object):
-    HGT_SIZE = (3601 * 3601 * 2)
-    
+class ElevationDownloader(object):
     def __init__(self, outpath='.'):
         self.outpath = outpath
 
@@ -27,6 +25,63 @@ def download_bboxes(self, bboxes):
         for name, bbox in bboxes.items():
             self.download_bbox(bbox)
     
+    def download_bbox(self, bbox, bucket='elevation-tiles-prod', prefix='geotiff'):
+        tiles = self.get_bbox_tiles(bbox)
+        found = set()
+        download = set()
+        for z,x,y in tiles:
+            od = self.tilepath(z, x, y)
+            op = os.path.join(self.outpath, *od)
+            if self.tile_exists(op):
+                pass
+                # found.add((x,y))
+            else:
+                download.add((x,y))
+        log.info("found %s tiles; %s to download"%(len(found), len(download)))
+        for x,y in sorted(download):
+            self.download_tile(bucket, prefix, z, x, y)
+
+    def tile_exists(self, op):
+        if os.path.exists(op):
+            return True
+
+    def download_tile(self, bucket, prefix, z, x, y, suffix=''):
+        od = self.tilepath(z, x, y)
+        op = os.path.join(self.outpath, *od)
+        makedirs(os.path.join(self.outpath, *od[:-1]))
+        if prefix:
+            od = [prefix]+od
+        url = 'http://s3.amazonaws.com/%s/%s%s'%(bucket, '/'.join(od), suffix)
+        log.info("downloading %s to %s"%(url, op))
+        download.download(url, op)
+        
+    def tilepath(self, z, x, y):
+        raise NotImplementedError
+
+    def get_bbox_tiles(self, bbox):
+        raise NotImplementedError
+
+class ElevationGeotiffDownloader(ElevationDownloader):
+    def get_bbox_tiles(self, bbox):
+        left, bottom, right, top = validate_bbox(bbox)
+        print "left, right, bottom, top:", left, right, bottom, top
+        zoom = 12
+        size = 2**zoom
+        xt = lambda x:int((x + 180.0) / 360.0 * size)
+        yt = lambda y:int((1.0 - math.log(math.tan(math.radians(y)) + (1 / math.cos(math.radians(y)))) / math.pi) / 2.0 * size)
+        tiles = []
+        for x in range(xt(left), xt(right)+1):
+            for y in range(yt(top), yt(bottom)+1):
+                print x, y
+                tiles.append([zoom, x, y])
+        return tiles
+
+    def tilepath(self, z, x, y):
+        return map(str, [z, x, str(y)+'.tif'])
+
+class ElevationHGTDownloader(ElevationDownloader):
+    HGT_SIZE = (3601 * 3601 * 2)
+    
     def get_bbox_tiles(self, bbox):
         left, bottom, right, top = validate_bbox(bbox)
         min_x = int(math.floor(left))
@@ -37,39 +92,15 @@ def get_bbox_tiles(self, bbox):
         tiles = set()
         for x in range(min_x, max_x):
             for y in range(min_y, max_y):
-                tiles.add((x,y))
+                tiles.add((0, x, y))
         return tiles
     
-    def download_bbox(self, bbox, bucket='elevation-tiles-prod', prefix='skadi'):
-        tiles = self.get_bbox_tiles(bbox)
-        found = set()
-        download = set()
-        for x,y in tiles:
-            od, key = self.hgtpath(x, y)
-            op = os.path.join(self.outpath, od, key)
-            if os.path.exists(op) and os.stat(op).st_size == self.HGT_SIZE:
-                found.add((x,y))
-            else:
-                download.add((x,y))
-        log.info("found %s tiles; %s to download"%(len(found), len(download)))
-        if len(download) > 100:
-            log.warning("  warning: downloading %s tiles will take an additional %0.2f GiB disk space"%(
-                len(download),
-                (len(download) * self.HGT_SIZE) / (1024.0**3)
-            ))
-        for x,y in sorted(download):
-            self.download_hgt(bucket, prefix, x, y)
-    
-    def hgtpath(self, x, y):
+    def download_tile(self, *args, **kwargs):
+        kwargs['prefix'] = 'skadi'
+        kwargs['suffix'] = '.gz'
+        super(*args, **kwargs)
+
+    def tilepath(self, z, x, y):
         ns = lambda i:'S%02d'%abs(i) if i < 0 else 'N%02d'%abs(i)
         ew = lambda i:'W%03d'%abs(i) if i < 0 else 'E%03d'%abs(i)
         return ns(y), '%s%s.hgt'%(ns(y), ew(x))
-
-    def download_hgt(self, bucket, prefix, x, y):
-        od, key = self.hgtpath(x, y)
-        op = os.path.join(self.outpath, od, key)
-        makedirs(os.path.join(self.outpath, od))
-        url = 'http://s3.amazonaws.com/%s/%s/%s/%s.gz'%(bucket, prefix, od, key)
-        log.info("downloading %s to %s"%(url, op))
-        download.download_gzip(url, op)
-        

From ff9a34c644d4df45aaccff5202b555751a9c84e1 Mon Sep 17 00:00:00 2001
From: Ian Rees <ian@ianrees.net>
Date: Wed, 28 Nov 2018 15:58:15 -0800
Subject: [PATCH 02/13] Options

---
 planetutils/elevation_tile_download.py   | 12 +++++++++++-
 planetutils/elevation_tile_downloader.py |  2 +-
 2 files changed, 12 insertions(+), 2 deletions(-)

diff --git a/planetutils/elevation_tile_download.py b/planetutils/elevation_tile_download.py
index e779ca3..aa2281b 100644
--- a/planetutils/elevation_tile_download.py
+++ b/planetutils/elevation_tile_download.py
@@ -12,12 +12,22 @@ def main():
     parser.add_argument('--csv', help='Path to CSV file with bounding box definitions.')
     parser.add_argument('--bbox', help='Bounding box for extract file. Format for coordinates: left,bottom,right,top')
     parser.add_argument('--verbose', help="Verbose output", action='store_true')
+    parser.add_argument('--format', help='Download format', default='geotiff')
+    parser.add_argument('--merge', help='Merge GeoTIFF tiles into output file')
+    parser.add_argument('--resample', help='Resample 16 bit GeoTIFF tiles into 8 bit with given min,max range')
     args = parser.parse_args()
 
     if args.verbose:
         log.set_verbose()
 
-    p = ElevationGeotiffDownloader(args.outpath)
+    if args.format == 'geotiff':
+        p = ElevationGeotiffDownloader(args.outpath)
+    elif args.format == 'skadi':
+        p = ElevationSkadiDownloader(args.outpath)
+    else:
+        print "Unknown format: %s"%args.format
+        sys.exit(1)
+        
     if args.csv:
         p.download_bboxes(load_bboxes_csv(args.csv))
     elif args.bbox:
diff --git a/planetutils/elevation_tile_downloader.py b/planetutils/elevation_tile_downloader.py
index 60c0aa3..e003726 100644
--- a/planetutils/elevation_tile_downloader.py
+++ b/planetutils/elevation_tile_downloader.py
@@ -79,7 +79,7 @@ def get_bbox_tiles(self, bbox):
     def tilepath(self, z, x, y):
         return map(str, [z, x, str(y)+'.tif'])
 
-class ElevationHGTDownloader(ElevationDownloader):
+class ElevationSkadiDownloader(ElevationDownloader):
     HGT_SIZE = (3601 * 3601 * 2)
     
     def get_bbox_tiles(self, bbox):

From 8135415c89adf2c9fe67d690a9cdbc81b90a435f Mon Sep 17 00:00:00 2001
From: Ian Rees <ian@ianrees.net>
Date: Wed, 28 Nov 2018 18:22:41 -0800
Subject: [PATCH 03/13] Fix tests

---
 planetutils/elevation_tile_download.py   | 11 +++---
 planetutils/elevation_tile_downloader.py | 44 ++++++++++++++----------
 tests/test_elevation_tile_downloader.py  | 25 +++++++-------
 3 files changed, 45 insertions(+), 35 deletions(-)

diff --git a/planetutils/elevation_tile_download.py b/planetutils/elevation_tile_download.py
index aa2281b..70caffa 100644
--- a/planetutils/elevation_tile_download.py
+++ b/planetutils/elevation_tile_download.py
@@ -4,7 +4,7 @@
 
 from . import log
 from .bbox import load_bboxes_csv, bbox_string
-from .elevation_tile_downloader import ElevationGeotiffDownloader
+from .elevation_tile_downloader import ElevationGeotiffDownloader, ElevationSkadiDownloader
 
 def main():
     parser = argparse.ArgumentParser()
@@ -13,21 +13,21 @@ def main():
     parser.add_argument('--bbox', help='Bounding box for extract file. Format for coordinates: left,bottom,right,top')
     parser.add_argument('--verbose', help="Verbose output", action='store_true')
     parser.add_argument('--format', help='Download format', default='geotiff')
-    parser.add_argument('--merge', help='Merge GeoTIFF tiles into output file')
-    parser.add_argument('--resample', help='Resample 16 bit GeoTIFF tiles into 8 bit with given min,max range')
+    parser.add_argument('--zoom', help='Zoom level', default=0, type=int)
+
     args = parser.parse_args()
 
     if args.verbose:
         log.set_verbose()
 
     if args.format == 'geotiff':
-        p = ElevationGeotiffDownloader(args.outpath)
+        p = ElevationGeotiffDownloader(args.outpath, zoom=args.zoom)
     elif args.format == 'skadi':
         p = ElevationSkadiDownloader(args.outpath)
     else:
         print "Unknown format: %s"%args.format
         sys.exit(1)
-        
+
     if args.csv:
         p.download_bboxes(load_bboxes_csv(args.csv))
     elif args.bbox:
@@ -35,5 +35,6 @@ def main():
     else:
         p.download_planet()
 
+
 if __name__ == '__main__':
     main()
diff --git a/planetutils/elevation_tile_downloader.py b/planetutils/elevation_tile_downloader.py
index e003726..b12de01 100644
--- a/planetutils/elevation_tile_downloader.py
+++ b/planetutils/elevation_tile_downloader.py
@@ -30,11 +30,10 @@ def download_bbox(self, bbox, bucket='elevation-tiles-prod', prefix='geotiff'):
         found = set()
         download = set()
         for z,x,y in tiles:
-            od = self.tilepath(z, x, y)
+            od = self.tile_path(z, x, y)
             op = os.path.join(self.outpath, *od)
             if self.tile_exists(op):
-                pass
-                # found.add((x,y))
+                found.add((x,y))
             else:
                 download.add((x,y))
         log.info("found %s tiles; %s to download"%(len(found), len(download)))
@@ -46,37 +45,41 @@ def tile_exists(self, op):
             return True
 
     def download_tile(self, bucket, prefix, z, x, y, suffix=''):
-        od = self.tilepath(z, x, y)
+        od = self.tile_path(z, x, y)
         op = os.path.join(self.outpath, *od)
         makedirs(os.path.join(self.outpath, *od[:-1]))
         if prefix:
             od = [prefix]+od
         url = 'http://s3.amazonaws.com/%s/%s%s'%(bucket, '/'.join(od), suffix)
         log.info("downloading %s to %s"%(url, op))
-        download.download(url, op)
+        self._download(url, op)
         
-    def tilepath(self, z, x, y):
+    def tile_path(self, z, x, y):
         raise NotImplementedError
 
     def get_bbox_tiles(self, bbox):
         raise NotImplementedError
 
+    def _download(url, op):
+        download.download(url, op)
+
 class ElevationGeotiffDownloader(ElevationDownloader):
+    def __init__(self, *args, **kwargs):
+        self.zoom = kwargs.pop('zoom', 0)
+        super(ElevationGeotiffDownloader, self).__init__(*args, **kwargs)
+
     def get_bbox_tiles(self, bbox):
         left, bottom, right, top = validate_bbox(bbox)
-        print "left, right, bottom, top:", left, right, bottom, top
-        zoom = 12
-        size = 2**zoom
+        size = 2**self.zoom
         xt = lambda x:int((x + 180.0) / 360.0 * size)
         yt = lambda y:int((1.0 - math.log(math.tan(math.radians(y)) + (1 / math.cos(math.radians(y)))) / math.pi) / 2.0 * size)
         tiles = []
         for x in range(xt(left), xt(right)+1):
             for y in range(yt(top), yt(bottom)+1):
-                print x, y
-                tiles.append([zoom, x, y])
+                tiles.append([self.zoom, x, y])
         return tiles
 
-    def tilepath(self, z, x, y):
+    def tile_path(self, z, x, y):
         return map(str, [z, x, str(y)+'.tif'])
 
 class ElevationSkadiDownloader(ElevationDownloader):
@@ -95,12 +98,17 @@ def get_bbox_tiles(self, bbox):
                 tiles.add((0, x, y))
         return tiles
     
-    def download_tile(self, *args, **kwargs):
-        kwargs['prefix'] = 'skadi'
-        kwargs['suffix'] = '.gz'
-        super(*args, **kwargs)
+    def tile_exists(self, op):
+        if os.path.exists(op) and os.stat(op).st_size == self.HGT_SIZE:	
+            return True
+
+    def download_tile(self, bucket, prefix, z, x, y, suffix=''):
+        super(ElevationSkadiDownloader, self).download_tile(bucket, 'skadi', z, x, y, suffix='.gz')
 
-    def tilepath(self, z, x, y):
+    def tile_path(self, z, x, y):
         ns = lambda i:'S%02d'%abs(i) if i < 0 else 'N%02d'%abs(i)
         ew = lambda i:'W%03d'%abs(i) if i < 0 else 'E%03d'%abs(i)
-        return ns(y), '%s%s.hgt'%(ns(y), ew(x))
+        return [ns(y), '%s%s.hgt'%(ns(y), ew(x))]
+
+    def _download(self, url, op):
+        download.download_gzip(url, op)
diff --git a/tests/test_elevation_tile_downloader.py b/tests/test_elevation_tile_downloader.py
index 1e8e45d..a18490a 100644
--- a/tests/test_elevation_tile_downloader.py
+++ b/tests/test_elevation_tile_downloader.py
@@ -4,7 +4,7 @@
 import types
 import unittest
 
-from planetutils.elevation_tile_downloader import ElevationTileDownloader
+from planetutils.elevation_tile_downloader import ElevationSkadiDownloader, ElevationGeotiffDownloader
 
 CA = [-126.386719,32.157012,-113.532715,42.244785]
 
@@ -13,14 +13,14 @@ def test_download_bboxes(self):
         pass
         
     def test_hgtpath(self):
-        e = ElevationTileDownloader('.')
+        e = ElevationSkadiDownloader('.')
         expect = ('N122', 'N122E037.hgt')
-        hgtpath = e.hgtpath(37, 122)
+        hgtpath = e.tile_path(0, 37, 122)
         self.assertEqual(hgtpath[0], expect[0])
         self.assertEqual(hgtpath[1], expect[1])
     
     def test_get_bbox_tiles(self):
-        e = ElevationTileDownloader('.')
+        e = ElevationSkadiDownloader('.')
         tiles = e.get_bbox_tiles(CA)
         self.assertEqual(len(tiles), 154)
         tiles = e.get_bbox_tiles([-180,-90,180,90])
@@ -28,31 +28,32 @@ def test_get_bbox_tiles(self):
     
     def download_bbox(self, e, method, args, expect):
         COUNT = []
-        def c(self, bucket, prefix, x, y):
-            COUNT.append([x,y])
-        e.download_hgt = types.MethodType(c, ElevationTileDownloader)
+        # def c(self, url, op):
+        def c(self, bucket, prefix, z, x, y, suffix=''):
+            COUNT.append([x, y])
+        e.download_tile = types.MethodType(c, ElevationSkadiDownloader)
         method(*args)
         self.assertEqual(len(COUNT), expect)
     
     def test_download_planet(self):
-        e = ElevationTileDownloader('.')
+        e = ElevationSkadiDownloader('.')
         self.download_bbox(e, e.download_planet, [], 64800)
     
     def test_download_bbox(self):
-        e = ElevationTileDownloader('.')
+        e = ElevationSkadiDownloader('.')
         self.download_bbox(e, e.download_bbox, [CA], 154)
     
     def test_download_bbox_found(self):
         d = tempfile.mkdtemp()
-        e = ElevationTileDownloader(d)
+        e = ElevationSkadiDownloader(d)
         # correct size
-        path = e.hgtpath(-119, 37)
+        path = e.tile_path(0, -119, 37)
         os.makedirs(os.path.join(d, path[0]))
         dp1 = os.path.join(d, *path)
         with open(dp1, 'w') as f:
             f.write('0'*e.HGT_SIZE)
         # incorrect size
-        path = e.hgtpath(-119, 36)
+        path = e.tile_path(0, -119, 36)
         os.makedirs(os.path.join(d, path[0]))
         dp2 = os.path.join(d, *path)
         with open(dp2, 'w') as f:

From 617cb2d2a5450272f00983f6ed6b6a3151c33daa Mon Sep 17 00:00:00 2001
From: Ian Rees <ian@ianrees.net>
Date: Wed, 28 Nov 2018 18:49:40 -0800
Subject: [PATCH 04/13] More tests

---
 planetutils/elevation_tile_downloader.py |  9 ++++++++-
 tests/test_elevation_tile_downloader.py  | 18 +++++++++++++++++-
 2 files changed, 25 insertions(+), 2 deletions(-)

diff --git a/planetutils/elevation_tile_downloader.py b/planetutils/elevation_tile_downloader.py
index b12de01..e34189c 100644
--- a/planetutils/elevation_tile_downloader.py
+++ b/planetutils/elevation_tile_downloader.py
@@ -70,10 +70,17 @@ def __init__(self, *args, **kwargs):
 
     def get_bbox_tiles(self, bbox):
         left, bottom, right, top = validate_bbox(bbox)
+        ybound = 85.0511
+        if bottom <= -ybound:
+            bottom = -ybound
+        if top > ybound:
+            top = ybound
+        if right >= 180:
+            right = 179.999
         size = 2**self.zoom
         xt = lambda x:int((x + 180.0) / 360.0 * size)
         yt = lambda y:int((1.0 - math.log(math.tan(math.radians(y)) + (1 / math.cos(math.radians(y)))) / math.pi) / 2.0 * size)
-        tiles = []
+        tiles = []    
         for x in range(xt(left), xt(right)+1):
             for y in range(yt(top), yt(bottom)+1):
                 tiles.append([self.zoom, x, y])
diff --git a/tests/test_elevation_tile_downloader.py b/tests/test_elevation_tile_downloader.py
index a18490a..29b6085 100644
--- a/tests/test_elevation_tile_downloader.py
+++ b/tests/test_elevation_tile_downloader.py
@@ -8,7 +8,23 @@
 
 CA = [-126.386719,32.157012,-113.532715,42.244785]
 
-class TestElevationTileDownloader(unittest.TestCase):
+class TestGeotiffDownloader(unittest.TestCase):
+    def test_tile_path(self):
+        e = ElevationGeotiffDownloader('.')
+        expect = ('0', '37', '122.tif')
+        tile_path = e.tile_path(0, 37, 122)
+        self.assertEqual(tile_path[0], expect[0])
+        self.assertEqual(tile_path[1], expect[1])
+        self.assertEqual(tile_path[2], expect[2])
+
+    def test_get_bbox_tiles(self):
+        e = ElevationGeotiffDownloader('.', zoom=8)
+        tiles = e.get_bbox_tiles(CA)
+        self.assertEqual(len(tiles), 100)
+        tiles = e.get_bbox_tiles([-180,-90,180,90])
+        self.assertEqual(len(tiles), 2**16)
+
+class TestElevationSkadiDownloader(unittest.TestCase):
     def test_download_bboxes(self):
         pass
         

From 7e70a2809bc0265c3522f17b9f2d942c7a6809e1 Mon Sep 17 00:00:00 2001
From: Ian Rees <ian@ianrees.net>
Date: Wed, 28 Nov 2018 19:36:32 -0800
Subject: [PATCH 05/13] elevation_tile_merge

---
 planetutils/elevation_tile_downloader.py |  2 +-
 planetutils/elevation_tile_merge.py      | 61 ++++++++++++++++++++++++
 setup.py                                 |  1 +
 3 files changed, 63 insertions(+), 1 deletion(-)
 create mode 100644 planetutils/elevation_tile_merge.py

diff --git a/planetutils/elevation_tile_downloader.py b/planetutils/elevation_tile_downloader.py
index e34189c..1e70d45 100644
--- a/planetutils/elevation_tile_downloader.py
+++ b/planetutils/elevation_tile_downloader.py
@@ -60,7 +60,7 @@ def tile_path(self, z, x, y):
     def get_bbox_tiles(self, bbox):
         raise NotImplementedError
 
-    def _download(url, op):
+    def _download(self, url, op):
         download.download(url, op)
 
 class ElevationGeotiffDownloader(ElevationDownloader):
diff --git a/planetutils/elevation_tile_merge.py b/planetutils/elevation_tile_merge.py
new file mode 100644
index 0000000..5a6bbfe
--- /dev/null
+++ b/planetutils/elevation_tile_merge.py
@@ -0,0 +1,61 @@
+#!/usr/bin/env python
+from __future__ import absolute_import, unicode_literals
+import argparse
+import sys
+import fnmatch
+import os
+import subprocess
+import tempfile
+
+from . import log
+
+def main():
+    parser = argparse.ArgumentParser()
+    parser.add_argument('--scale', help="Resample to 8 bit with (min,max) range")
+    parser.add_argument('outpath', help='Output filename')
+    parser.add_argument('inpath', help='Input directory')
+    args = parser.parse_args()
+
+    outpath = args.outpath
+    tmppath = args.outpath
+
+    if args.scale and len(args.scale.split(',')) != 2:
+        print "Must provide min, max values"
+        sys.exit(1)
+    elif args.scale:
+        # Output to tmp file
+        _, tmppath = tempfile.mkstemp(suffix='.tif')
+
+    matches = []
+    for root, dirnames, filenames in os.walk(args.inpath):
+        for filename in fnmatch.filter(filenames, '*.tif'):
+            matches.append(os.path.join(root, filename))
+
+    if len(matches) == 0:
+        print "No input files"
+        sys.exit(0)
+
+    print "Found %s files:"%len(matches)
+    for i in matches:
+        print "\t%s"%(i)
+
+    # gdal_merge.py -init 0 -o out.tif
+    print "Merging... %s"%(tmppath)
+    cmd = ['gdal_merge.py', '-init', '0', '-o', tmppath]
+    cmd += matches
+    p = subprocess.check_call(cmd)
+
+    # gdal_translate -of GTiff -ot Byte -scale 0 255 0 255 out.tif out8.tif
+    if args.scale:
+        print "Scaling: %s -> %s"%(tmppath, outpath)
+        a = args.scale.split(",")
+        cmd = ['gdal_translate', '-of', 'GTiff', '-ot', 'Byte', '-scale', a[0], a[1], '0', '255', tmppath, outpath]
+        subprocess.check_call(cmd)
+        # cleanup
+        try: os.unlink('%s.aux.xml'%outpath)
+        except: pass
+        try: os.unlink(tmppath)
+        except: pass
+
+if __name__ == '__main__':
+    main()
diff --git a/setup.py b/setup.py
index 5a759c9..ca2a498 100644
--- a/setup.py
+++ b/setup.py
@@ -27,6 +27,7 @@
             'osm_planet_get_timestamp=planetutils.osm_planet_get_timestamp:main',
             'osm_extract_download=planetutils.osm_extract_download:main',
             'elevation_tile_download=planetutils.elevation_tile_download:main',
+            'elevation_tile_merge=planetutils.elevation_tile_merge:main',
             'valhalla_tilepack_download=planetutils.tilepack_download:main',
             'valhalla_tilepack_list=planetutils.tilepack_list:main'
         ],

From d337c01353d53c9e2aa0e1d1b49cfcd1f07c66e9 Mon Sep 17 00:00:00 2001
From: Drew Dara-Abrams <drew@interline.io>
Date: Fri, 30 Nov 2018 14:23:30 -0800
Subject: [PATCH 06/13] CircleCI: try running tests within container

---
 .circleci/config.yml | 15 ++++++---------
 1 file changed, 6 insertions(+), 9 deletions(-)

diff --git a/.circleci/config.yml b/.circleci/config.yml
index c8d873b..e59b2c5 100644
--- a/.circleci/config.yml
+++ b/.circleci/config.yml
@@ -1,16 +1,13 @@
 version: 2
 jobs:
   build:
-    docker:
-      - image: circleci/python:3.6.6
+    machine: true
     steps:
       - checkout
+      - setup_remote_docker
       - run:
-          name: Installing OSM dependencies
-          command: sudo apt-get install osmosis osmctools osmium-tool libboost-python-dev libexpat1-dev zlib1g-dev libbz2-dev
+          name: Build container
+          command: docker build -t interline-io/planetutils:$CIRCLE_BRANCH .
       - run:
-          name: Installing Python package and dependencies
-          command: pip install --user .
-      - run:
-          name: Running tests
-          command: python setup.py test
+          name: Run tests in container
+          command: docker run -it interline-io/planetutils:$CIRCLE_BRANCH python setup.py test

From 1c3dcf67ff9b3b2c0298b43d62baa084c11e74b6 Mon Sep 17 00:00:00 2001
From: Drew Dara-Abrams <drew@interline.io>
Date: Fri, 30 Nov 2018 14:25:07 -0800
Subject: [PATCH 07/13] "You can use Docker without a special step"

---
 .circleci/config.yml | 1 -
 1 file changed, 1 deletion(-)

diff --git a/.circleci/config.yml b/.circleci/config.yml
index e59b2c5..7fdb8d3 100644
--- a/.circleci/config.yml
+++ b/.circleci/config.yml
@@ -4,7 +4,6 @@ jobs:
     machine: true
     steps:
       - checkout
-      - setup_remote_docker
       - run:
           name: Build container
           command: docker build -t interline-io/planetutils:$CIRCLE_BRANCH .

From 8e1bfb89aa7bf9353aaad887b99c665ca9fdd9ef Mon Sep 17 00:00:00 2001
From: Drew Dara-Abrams <drew@interline.io>
Date: Fri, 30 Nov 2018 14:29:40 -0800
Subject: [PATCH 08/13] building container will run tests

---
 .circleci/config.yml | 5 +----
 1 file changed, 1 insertion(+), 4 deletions(-)

diff --git a/.circleci/config.yml b/.circleci/config.yml
index 7fdb8d3..ca349a4 100644
--- a/.circleci/config.yml
+++ b/.circleci/config.yml
@@ -5,8 +5,5 @@ jobs:
     steps:
       - checkout
       - run:
-          name: Build container
+          name: Build container (which includes running tests)
           command: docker build -t interline-io/planetutils:$CIRCLE_BRANCH .
-      - run:
-          name: Run tests in container
-          command: docker run -it interline-io/planetutils:$CIRCLE_BRANCH python setup.py test

From 172ec507ff7c52a04a1825bbe619ac161f81b126 Mon Sep 17 00:00:00 2001
From: Drew Dara-Abrams <drew@interline.io>
Date: Fri, 30 Nov 2018 14:42:57 -0800
Subject: [PATCH 09/13] CircleCI: you should fail on this!

---
 tests/test_bbox.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/tests/test_bbox.py b/tests/test_bbox.py
index 1212e5d..56bcac3 100644
--- a/tests/test_bbox.py
+++ b/tests/test_bbox.py
@@ -22,9 +22,9 @@ def test_bounds(self):
         self.assertRaises(AssertionError, bbox.validate_bbox, (0,0,-180,0))
         self.assertRaises(AssertionError, bbox.validate_bbox, (0,90,0,0))
         self.assertRaises(AssertionError, bbox.validate_bbox, (0,0,0,-90))
-    
+
     def test_returns_array(self):
-        self.assertEqual(bbox.validate_bbox([1,2,3,4]), [1.0, 2.0, 3.0, 4.0])
+        self.assertEqual(bbox.validate_bbox([1,2,3]), [1.0, 2.0, 3.0, 4.0])
 
 class TestLoadBboxesCsv(unittest.TestCase):
     def test_load(self):
@@ -39,7 +39,7 @@ def test_load(self):
 class TestBboxString(unittest.TestCase):
     def test_returns_array(self):
         self.assertEqual(bbox.bbox_string('1.0,2.0,3.0,4.0'), [1.0,2.0,3.0,4.0])
-    
+
     def test_validates(self):
         self.assertRaises(AssertionError, bbox.bbox_string, ('10,-10,20,-20'))
 

From 660c70d3e2d8985c34e24aaf04761c17acc88dfb Mon Sep 17 00:00:00 2001
From: Drew Dara-Abrams <drew@interline.io>
Date: Fri, 30 Nov 2018 14:46:55 -0800
Subject: [PATCH 10/13] CircleCI: now you should pass!

---
 tests/test_bbox.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/test_bbox.py b/tests/test_bbox.py
index 56bcac3..4791034 100644
--- a/tests/test_bbox.py
+++ b/tests/test_bbox.py
@@ -24,7 +24,7 @@ def test_bounds(self):
         self.assertRaises(AssertionError, bbox.validate_bbox, (0,0,0,-90))
 
     def test_returns_array(self):
-        self.assertEqual(bbox.validate_bbox([1,2,3]), [1.0, 2.0, 3.0, 4.0])
+        self.assertEqual(bbox.validate_bbox([1,2,3,4]), [1.0, 2.0, 3.0, 4.0])
 
 class TestLoadBboxesCsv(unittest.TestCase):
     def test_load(self):

From 966053e4477eb8c3e3428ad8e8c90a0344d58bfb Mon Sep 17 00:00:00 2001
From: Drew Dara-Abrams <drew@interline.io>
Date: Fri, 30 Nov 2018 15:03:49 -0800
Subject: [PATCH 11/13] add GDAL dependency

---
 Dockerfile | 1 +
 README.md  | 1 +
 2 files changed, 2 insertions(+)

diff --git a/Dockerfile b/Dockerfile
index 02afb5a..e4c4ee1 100644
--- a/Dockerfile
+++ b/Dockerfile
@@ -11,6 +11,7 @@ RUN apt-get update -y && apt-get install \
       osmctools \
       osmium-tool \
       pyosmium \
+      gdal \
       awscli \
       software-properties-common \
       -y
diff --git a/README.md b/README.md
index 934ae98..a36947e 100644
--- a/README.md
+++ b/README.md
@@ -85,6 +85,7 @@ If you want to install and use the Python package directly, you'll need to provi
 - [OSM C tools](https://gitlab.com/osm-c-tools/osmctools/)
 - [Osmium Tool](https://osmcode.org/osmium-tool/)
 - [PyOsmium](https://osmcode.org/pyosmium/)
+- [GDAL](https://www.gdal.org/)
 
 Then clone this repo, run the tests, and install the Python package:
 

From e70e052c4826c2f94e72a9c776921a034aee5e0f Mon Sep 17 00:00:00 2001
From: Drew Dara-Abrams <drew@interline.io>
Date: Fri, 30 Nov 2018 15:07:41 -0800
Subject: [PATCH 12/13] readme: additions

---
 README.md | 32 +++++++++++++++++++++++++++++---
 1 file changed, 29 insertions(+), 3 deletions(-)

diff --git a/README.md b/README.md
index a36947e..00bbbb9 100644
--- a/README.md
+++ b/README.md
@@ -22,6 +22,7 @@
   * [osm_extract_download](#osm_extract_download)
   * [osm_planet_get_timestamp](#osm_planet_get_timestamp)
   * [elevation_tile_download](#elevation_tile_download)
+  * [elevation_tile_merge](#elevation_tile_merge)
   * [valhalla_tilepack_list](#valhalla_tilepack_list)
   * [valhalla_tilepack_download](#valhalla_tilepack_download)
 - [Specifying bounding boxes](#specifying-bounding-boxes)
@@ -40,6 +41,7 @@ Python-based scripts and a Docker container to work with planet-scale geographic
 - cut your copy of the OSM planet into named bounding boxes
 - download [OSM Extracts from Interline](https://www.interline.io/osm/extracts/) for popular cities and regions
 - download [Mapzen Terrain Tiles from AWS](https://aws.amazon.com/public-datasets/terrain/) for the planet or your bounding boxes
+- merge and resample Terrain Tiles
 - download [Valhalla Tilepacks from Interline](https://www.interline.io/valhalla/tilepacks) for the planet (subscription required)
 
 PlanetUtils is packaged for use as a:
@@ -174,13 +176,17 @@ osm_planet_get_timestamp planet-latest.osm.pbf
 
 Download elevation tiles from the [Terrain Tiles in the AWS Public Datasets program](https://aws.amazon.com/public-datasets/terrain/). Download for the entire planet, only tiles within a single bounding box, or within multiple bounding boxes.
 
-To download the entire planet of tiles (__which will require about 1.6Tb of space!__):
+Elevation tiles are available in [a variety of formats](https://mapzen.com/documentation/terrain-tiles/formats/). This command supports the download of:
+- GeoTIFF (default): extension `.tif` in Web Mercator projection, 512x512 tiles
+- Skadi: extension `.hgt` in unprojected latlng, 1°x1° tiles
+
+To download the entire planet in Skadi tiles (__which will require about 1.6Tb of space!__):
 
 ```sh
-elevation_tile_download --outpath=data/elevation
+elevation_tile_download --format=skadi --outpath=data/elevation
 ```
 
-To download tiles to cover a single bounding box:
+To download GeoTIFF tiles to cover a single bounding box:
 
 ```sh
 elevation_tile_download --outpath=data/elevation --bbox=-122.737,37.449,-122.011,37.955
@@ -198,6 +204,26 @@ For complete help on command-line arguments:
 elevation_tile_download -h
 ```
 
+### elevation_tile_merge
+
+After downloading elevation tiles using the `elevation_tile_download` command, use this command to merge together multiple tiles. You can optionally resample elevation values as part of the merge process.
+
+This command only operates on GeoTIFF format elevation tiles.
+
+Warnings: merging lots of tiles can be resource intensive!
+
+To merge a directory of GeoTIFF files:
+
+```sh
+elevation_tile_merge geo_tiff_tiles/ single_tile.tif
+```
+
+For complete help on command-line arguments:
+
+```sh
+elevation_tile_merge -h
+```
+
 ### valhalla_tilepack_list
 
 Use [Valhalla Tilepacks from Interline](https://www.interline.io/valhalla/tilepacks/) to power your own instances of the [Valhalla routing engine](https://www.interline.io/valhalla/). Anyone can list available planet tilepacks. A subscription and an API key are required to [download tilepacks](#valhalla_tilepack_download).

From 2e7947b0980ca6a80d492cc38caae8d0599d1e99 Mon Sep 17 00:00:00 2001
From: Drew Dara-Abrams <drew@interline.io>
Date: Fri, 30 Nov 2018 15:12:56 -0800
Subject: [PATCH 13/13] different GDAL pacakge on Ubuntu

---
 Dockerfile | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Dockerfile b/Dockerfile
index e4c4ee1..8fae76c 100644
--- a/Dockerfile
+++ b/Dockerfile
@@ -11,7 +11,7 @@ RUN apt-get update -y && apt-get install \
       osmctools \
       osmium-tool \
       pyosmium \
-      gdal \
+      python-gdal \
       awscli \
       software-properties-common \
       -y