Merge pull request #22 from interline-io/download-geotif

Additional elevation tile download options

.circleci/config.yml

@@ -1,16 +1,9 @@
 version: 2
 jobs:
   build:
-    docker:
-      - image: circleci/python:3.6.6
+    machine: true
     steps:
       - checkout
       - run:
-          name: Installing OSM dependencies
-          command: sudo apt-get install osmosis osmctools osmium-tool libboost-python-dev libexpat1-dev zlib1g-dev libbz2-dev
-      - run:
-          name: Installing Python package and dependencies
-          command: pip install --user .
-      - run:
-          name: Running tests
-          command: python setup.py test
+          name: Build container (which includes running tests)
+          command: docker build -t interline-io/planetutils:$CIRCLE_BRANCH .

Dockerfile

@@ -11,6 +11,7 @@ RUN apt-get update -y && apt-get install \
       osmctools \
       osmium-tool \
       pyosmium \
+      python-gdal \
       awscli \
       software-properties-common \
       -y

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:
@@ -85,6 +87,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:
 
@@ -173,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
@@ -197,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).

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)

planetutils/elevation_tile_download.py

@@ -4,26 +4,37 @@
 
 from . import log
 from .bbox import load_bboxes_csv, bbox_string
-from .elevation_tile_downloader import ElevationTileDownloader
+from .elevation_tile_downloader import ElevationGeotiffDownloader, ElevationSkadiDownloader
 
 def main():
     parser = argparse.ArgumentParser()
     parser.add_argument('--outpath', help='Output path for elevation tiles.', default='.')
     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('--zoom', help='Zoom level', default=0, type=int)
+
     args = parser.parse_args()
 
     if args.verbose:
         log.set_verbose()
 
-    p = ElevationTileDownloader(args.outpath)
+    if args.format == 'geotiff':
+        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:
         p.download_bbox(bbox_string(args.bbox))
     else:
         p.download_planet()
 
+
 if __name__ == '__main__':
     main()

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,73 @@ 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.tile_path(z, x, y)
+            op = os.path.join(self.outpath, *od)
+            if self.tile_exists(op):
+                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.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))
+        self._download(url, op)
+
+    def tile_path(self, z, x, y):
+        raise NotImplementedError
+
+    def get_bbox_tiles(self, bbox):
+        raise NotImplementedError
+
+    def _download(self, 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)
+        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 = []    
+        for x in range(xt(left), xt(right)+1):
+            for y in range(yt(top), yt(bottom)+1):
+                tiles.append([self.zoom, x, y])
+        return tiles
+
+    def tile_path(self, z, x, y):
+        return map(str, [z, x, str(y)+'.tif'])
+
+class ElevationSkadiDownloader(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 +102,20 @@ 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 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 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_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))
+    def _download(self, url, op):
         download.download_gzip(url, op)
-

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()

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'
         ],

tests/test_bbox.py

@@ -22,7 +22,7 @@ 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])
 
@@ -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'))