Adds support for processing raster data

This adds a `normalize` function that accepts a GeoTIFF file as input
and creates a new GeoTIFF which has been optimized for serving through
GeoServer. The following steps are taken:

* Paletted images are expanded to RGB.
* Images are internally tiled using a block size of 512.
* Overview images are added to the TIFF when larger than block size.

Support is provided through the `rasterio` package which itself depends
on GDAL. You'll need to make sure the header files for GDAL are
accessible in order to install it.

.travis.yml

@@ -5,7 +5,8 @@ services:
     - postgresql
 env:
   global:
-    POSTGIS_DB="postgresql://postgres@localhost/slingshot_test"
+    - POSTGIS_DB="postgresql://postgres@localhost/slingshot_test"
+    - GDAL_VERSION=2.2.0
 matrix:
   include:
     - env: TOX_ENV=py27-integration
@@ -15,9 +16,17 @@ matrix:
       env: TOX_ENV=flake8
     - python: 3.5
       env: TOX_ENV=coveralls
+cache:
+  directories:
+    - $HOME/gdal
+before_install:
+  - bash install-gdal.sh
+  - export LD_LIBRARY_PATH=$HOME/gdal/lib:$LD_LIBRARY_PATH
+  - export PATH=$HOME/gdal/bin:$PATH
 before_script:
-    - psql -U postgres -c "CREATE DATABASE slingshot_test;"
-    - psql -U postgres -d slingshot_test -c "CREATE EXTENSION postgis;"
+  - echo "Using gdal version $(gdal-config --version)"
+  - psql -U postgres -c "CREATE DATABASE slingshot_test;"
+  - psql -U postgres -d slingshot_test -c "CREATE EXTENSION postgis;"
 install:
   - pip install tox
 script:

Pipfile

@@ -7,11 +7,13 @@ arrow = "*"
 bagit = "*"
 click = "*"
 geomet = "*"
+numpy = "*"
 psycopg2 = "*"
 pyshp = "*"
 requests = "*"
 GeoAlchemy2 = "*"
 PlyPlus = "*"
+rasterio = "*"
 
 [dev-packages]
 bumpversion = "*"

install-gdal.sh

@@ -0,0 +1,10 @@
+#!/usr/bin/env sh
+set -e
+
+if [ ! -d "$HOME/gdal/lib" ]; then
+  wget http://download.osgeo.org/gdal/$GDAL_VERSION/gdal-$GDAL_VERSION.tar.xz
+  tar -xf gdal-$GDAL_VERSION.tar.xz
+  cd gdal-$GDAL_VERSION && ./configure --prefix=$HOME/gdal && make && make install
+else
+  echo "Using cached directory."
+fi

setup.py

@@ -34,12 +34,15 @@
         'click',
         'geoalchemy2',
         'geomet',
+        'numpy',
         'plyplus',
         'psycopg2',
         'pyshp',
+        'rasterio',
         'requests',
     ],
     setup_requires=[
+        'numpy', # this is required (but not specified) by rasterio
         'pytest-runner',
     ],
     tests_require=[

slingshot/raster.py

@@ -0,0 +1,94 @@
+import math
+
+import rasterio
+from rasterio.enums import ColorInterp, Resampling
+import numpy as np
+
+
+BLOCKSIZE = 512
+TIFF_DEFAULTS = {
+    'compress': 'JPEG',
+    'interleave': 'PIXEL',
+}
+
+
+def make_palette(cmap):
+    """Turn a color map dictionary into a numpy array."""
+    palette = np.zeros((len(cmap),), dtype=[('R', np.uint8), ('G', np.uint8),
+                                            ('B', np.uint8), ('A', np.uint8)])
+    for k, v in cmap.items():
+        palette[k] = v
+    return palette
+
+
+def windows(fp):
+    """Generate tuples of TIFF windows and blocks.
+
+    ``windows`` is a generator that produces tuples of ``rasterio``
+    windows and block generators. A block generator produces data arrays
+    for each band in the given window.
+
+    This can be used to iterate over chunks of large TIFFs without
+    reading the whole thing into memory at once.
+    """
+    try:
+        palette = make_palette(fp.colormap(1))
+    except:
+        palette = None
+    for _, window in fp.block_windows():
+        yield window, blocks(fp, window, palette)
+
+
+def blocks(fp, window, palette=None):
+    """Generate blocks of band data for a given window.
+
+    This produces numpy arrays of data for an image's bands in the given
+    window. Paletted images will be expanded to RGB bands.
+    """
+    cinterp = fp.colorinterp(1)
+    if cinterp == ColorInterp.palette:
+        block = fp.read(1, window=window)
+        expanded = np.take(palette, block)
+        for band in ('R', 'G', 'B'):
+            yield expanded[band]
+    elif cinterp == ColorInterp.gray:
+        yield fp.read(1, window=window)
+    else:
+        for band in range(1, 4):
+            yield fp.read(band, window=window)
+
+
+def overviews(size, blocksize):
+    """Compute optimal list of overview levels.
+
+    The size parameter should be max(width, height).
+    """
+    num_levels = int(math.ceil(math.log((size/blocksize), 2)))
+    return [2**y for y in range(1, num_levels + 1)]
+
+
+def normalize(filein, fileout):
+    """Create a normalized GeoTIFF file.
+
+    This will take the input GeoTIFF and create a GeoTIFF that is
+    optimized for serving through GeoServer. Compression is set to
+    JPEG, the TIFF is tiled using a block size of 512, and overviews
+    are added. Paletted TIFFs are expanded to 3-band images and RGB
+    images are converted to YCbCr color space.
+    """
+    with rasterio.open(filein, 'r') as fp_in:
+        kwargs = fp_in.profile
+        kwargs.update(TIFF_DEFAULTS)
+        if max(fp_in.width, fp_in.height) >= BLOCKSIZE:
+            kwargs.update({'tiled': True, 'blockxsize': BLOCKSIZE,
+                           'blockysize': BLOCKSIZE})
+        if fp_in.colorinterp(1) == ColorInterp.palette or fp_in.count == 3:
+            kwargs.update({'photometric': 'YCBCR', 'count': 3})
+        with rasterio.open(fileout, 'w', **kwargs) as fp_out:
+            for window, bands in windows(fp_in):
+                for band, block in enumerate(bands, start=1):
+                    fp_out.write(block, window=window, indexes=band)
+        with rasterio.open(fileout, 'r+') as fp:
+            factors = overviews(max(fp.width, fp.height), BLOCKSIZE)
+            fp.build_overviews(factors, Resampling.average)
+            fp.update_tags(ns='rio_overview', resampling='average')

tests/conftest.py

@@ -68,6 +68,21 @@ def fgdc():
     return _data_file('fixtures/bermuda/data/bermuda.xml')
 
 
+@pytest.fixture
+def rgb():
+    return _data_file('fixtures/rgb.tif')
+
+
+@pytest.fixture
+def grayscale():
+    return _data_file('fixtures/grayscale.tif')
+
+
+@pytest.fixture
+def paletted():
+    return _data_file('fixtures/paletted.tif')
+
+
 def _data_file(name):
     cur_dir = os.path.dirname(os.path.realpath(__file__))
     return os.path.join(cur_dir, name)

tests/test_raster.py

@@ -0,0 +1,37 @@
+from tempfile import NamedTemporaryFile
+
+import rasterio
+
+from slingshot.raster import normalize, BLOCKSIZE, overviews
+
+
+def test_normalize_tiles_layer(rgb):
+    with NamedTemporaryFile() as out:
+        normalize(rgb, out.name)
+        with rasterio.open(out.name, 'r') as fp:
+            assert fp.block_shapes[0] == (BLOCKSIZE, BLOCKSIZE)
+
+
+def test_normalize_turns_paletted_layer_into_rgb(paletted):
+    with NamedTemporaryFile() as out:
+        normalize(paletted, out.name)
+        with rasterio.open(out.name, 'r') as fp:
+            assert fp.count == 3
+
+
+def test_normalize_keeps_grayscale_as_grayscale(grayscale):
+    with NamedTemporaryFile() as out:
+        normalize(grayscale, out.name)
+        with rasterio.open(out.name, 'r') as fp:
+            assert fp.count == 1
+
+
+def test_overviews_generates_factors_of_overviews():
+    assert overviews(2048, 512) == [2, 4]
+
+
+def test_normalize_creates_overviews(rgb):
+    with NamedTemporaryFile() as out:
+        normalize(rgb, out.name)
+        with rasterio.open(out.name, 'r') as fp:
+            assert fp.overviews(1) == [2]

tox.ini

@@ -6,6 +6,7 @@ passenv = HOME POSTGIS_DB
 deps =
     pipenv
     python-dotenv
+    numpy
     py27: mock
     {coverage,coveralls}: pytest-cov
     coveralls: coveralls