Create a color relief by using a lookup table

color_relief.py

#!/usr/bin/env python3
from PIL import Image
from osgeo import gdal
gdal.UseExceptions()

import numpy as np
import gc

# ----------------------------------
# COLOR LOOKUP TABLE (gradient color) 256x256 8bit
# x-axis == greyscale value of hillshade (256 possibilities)
# y-axis == interpolated height value from DEM (256 possibilities)

# E +-----------------------------+
# L |                             |
# E |                             |
# V |                             |
# A |                             |
# T |                             |
# I |                             |
# O |                             |
# N |                             |
#   |                             |
# D |                             |
# E |                             |
# M +-----------------------------+
#   0  GREYSCALE VALUE HILLSHADE  256

# see: http://cartography.oregonstate.edu/pdf/2006_JennyHurni_SwissStyleShading.pdf
# see: Swiss-Style Colour Relief Shading Modulated by Elevation and by Exposure to Illumination
# ----------------------------------

# http://www.geeks3d.com/20100930/tutorial-first-steps-with-pil-python-imaging-library/
img = Image.open('lookup_table.png', 'r')
# http://code.activestate.com/recipes/577591-conversion-of-pil-image-and-numpy-array/
color_lookup = np.array(img)

# free memory
img = None
gc.collect()

# ----------------------------------
# DEM DATA ARRAY
# use it for looking up the y-axis (rows) in color_lookup table
# ----------------------------------

ds = gdal.Open('cubicspline_1x1/xeis_elev_resampled.flt')
#ds = gdal.Open('filtered/xeis_elev_filtered.flt')
band = ds.GetRasterBand(1)
# convert float32 to integers (for array index lookup)
color_lookup_row_idx = band.ReadAsArray().astype('int16')

# http://gis.stackexchange.com/a/146989/58989
alt_min, alt_max, alt_mean, alt_stddev = band.GetStatistics(False, True)

print(np.amax(color_lookup_row_idx), alt_max)
# http://stackoverflow.com/questions/6824122/mapping-a-numpy-array-in-place
# http://opencvpython.blogspot.co.at/2012/06/fast-array-manipulation-in-numpy.html
from numpy import rint
def ufunc_ed(array_):
    # clamp the height between the color lookup values of 0-255
    multiplied = array_ / alt_max * 255
    array_[:] = multiplied
    #array_[:] = rint(multiplied)
# clamp values
ufunc_ed(color_lookup_row_idx)

# free memory
ds = None
band = None
gc.collect()

# ----------------------------------
# HILLSHADE DATA ARRAY (grey values) 256 values, 8bit
# ----------------------------------

img = Image.open('cubicspline_1x1/xeis_hillshade_resampled.tif', 'r')
#img = Image.open('filtered/xeis_elev_hillshade.tif', 'r')
color_lookup_col_idx = np.array(img)

# checks
assert color_lookup_row_idx.shape == color_lookup_col_idx.shape
assert int(np.amax(color_lookup_row_idx)) <= 255
assert int(np.amax(color_lookup_col_idx)) <= 255

# looking up the color in color lookup table
# see: http://stackoverflow.com/a/35941469/1230358
result_data = color_lookup[color_lookup_col_idx, color_lookup_row_idx]
#result_data = color_lookup[color_lookup_row_idx, color_lookup_col_idx]

img = Image.fromarray(result_data, 'RGB')
img.save('color_relief.png')

#  brew install gdal --with-python --with-python3 --with-complete --with-libkml --with-mysql
#  python3 -m cProfile -s tottime colorramp.py 

lookup_table.png

lookuptable.py

#!/usr/bin/env python3
import numpy
from scipy import interpolate
from scipy.ndimage import interpolation

# The base idea is from a paper of Bernhard Jenny and Lorenz Hurni
# see: http://cartography.oregonstate.edu/pdf/2006_JennyHurni_SwissStyleShading.pdf

# The coding is mostly inspired by
# see: http://stackoverflow.com/questions/39489089/interpolating-elements-of-a-color-matrix-on-the-basis-of-some-given-reference-el
# see: http://stackoverflow.com/questions/39485178/two-dimensional-color-ramp-256x256-matrix-interpolated-from-4-corner-colors


# also see D3.js and color interpolation (LAB and HSC)
# see: http://bl.ocks.org/syntagmatic/5bbf30e8a658bcd5152b
# see: http://stackoverflow.com/a/37503031/1230358
# see: http://dsp.stackexchange.com/questions/13697/how-do-you-interpolate-between-points-in-an-image-2d-e-g-using-splines


DEBUG = True


if DEBUG:
    from matplotlib import pyplot
    from matplotlib import gridspec

    gs = gridspec.GridSpec(4, 3) # rows, cols
    
    font_size = 8
    pyplot.rc('font', size=font_size)        # controls default text sizes
    pyplot.rc('axes', titlesize=font_size)   # fontsize of the axes title
    pyplot.rc('axes', labelsize=font_size)   # fontsize of the x any y labels
    pyplot.rc('xtick', labelsize=font_size)  # fontsize of the tick labels
    pyplot.rc('ytick', labelsize=font_size)  # fontsize of the tick labels
    pyplot.rc('legend', fontsize=font_size)  # legend fontsize
    pyplot.rc('figure', titlesize=font_size)
    
    fig = pyplot.figure()
    
    subplot_id = 0
    
    def add_plot(matrix, title):
        global subplot_id

        ax = fig.add_subplot(gs[subplot_id])
        #ax = pyplot.subplot(pyplot_idx)
        ax.set_title(title)
        ax.imshow(numpy.array(matrix, dtype=numpy.uint8), interpolation='nearest')
    
        subplot_id += 1


nan = numpy.nan

# the matrix with the reference color elements filled with nan's
ref = numpy.empty([7, 7, 3]) * nan
ref[0][6] = (239,238,185)
ref[1][1] = (120,131,125)
ref[4][6] = (184,191,171)
#ref[6][2] = (150,168,158)
ref[6][2] = (255,100,100)
ref[6][5] = (166,180,166)


def interpolate_reference_points():
    """Linear interpolation with the given reference points"""
    points = numpy.where(numpy.isfinite(ref[:,:,:]))
    values = ref[points]
    
    shape = ref.shape
    grid_x, grid_y, grid_z = numpy.mgrid[0:shape[0], 0:shape[1], 0:shape[2]]

    # linear interpolation
    return interpolate.griddata(points, values, (grid_x, grid_y, grid_z), method='linear')
    
filled_grid = interpolate_reference_points()


if DEBUG:
    add_plot(ref, "Original reference points")
    add_plot(filled_grid, "Linear interpolation with the given reference points")


def calc_rgb(rp0, rp1, point_idx):
    """Extrapolate the RGB value of the new point in relation to the two reference points"""
    idx_0, r0, g0, b0 = rp0
    idx_1, r1, g1, b1 = rp1
    
    # modify for calculation
    point_idx += 1
    idx_0 += 1
    idx_1 += 1

    r = r0+(r1-r0)*(point_idx-idx_0)/(idx_1-idx_0)
    g = g0+(g1-g0)*(point_idx-idx_0)/(idx_1-idx_0)
    b = b0+(b1-b0)*(point_idx-idx_0)/(idx_1-idx_0)

    # constrain the values for rgb compatible values
    return numpy.clip((r,g,b), 0, 255)
    
    
def sort_axis_by_filled_values(matrix, axis=0):
    """Return the distinct rows or columns sorted descending by the sum of their allready calculated rgb values
    axis 0 == vertically (y-axis)
    axis 1 == horizontally (x-axis)
    """
    if axis == 0:
        filled, _ = numpy.where(numpy.isfinite(matrix[:,:,0]))
    else:
        _, filled = numpy.where(numpy.isfinite(matrix[:,:,0]))
        
    filled_indices = numpy.bincount(filled)

    # do not process allready fully filled x-axis elements (columns)
    max_elements = matrix.shape[1] - 1
    while max_elements >= 0:
        max_idx = numpy.where((filled_indices==max_elements))
        max_elements -= 1
        # nothing found, decrease maximum
        if max_idx[0].size == 0:
            continue
        # found a result, yield it
        yield max_idx[0]
        
        
def extrapolate_axis(matrix, axis=0, output_matrix=None):
    # create a the result matrix with the same shape as original matrix
    if output_matrix is None:
        # inline modifications to original matrix
        output_matrix = matrix
    else:
        # must have the same shape
        assert matrix.shape == output_matrix.shape
    
    for ids in sort_axis_by_filled_values(matrix, axis):
        for idx in ids:
            if axis == 0:
                sequence = matrix[idx] # rows
            else:
                sequence = matrix[:,idx] # cols
            
            filled_points = numpy.where(numpy.isfinite(sequence[:,0]))[0]
            
            if len(filled_points) <= 1:
                print("Zero or one item. Too less to extrapolate a RGB color gradient on this axis")
                return
        
            # get the first and the last element of the reference colors
            first, last = filled_points[[0,-1]]
    
            max_dim = filled_grid.shape[axis]
            t = numpy.linspace(0, max_dim-1, num=max_dim, dtype=numpy.uint8)
            nan_indices = numpy.delete(t, filled_points, 0)

            for nan_idx in nan_indices:
                # get the nan col_idx for this row
                rgb = calc_rgb(
                    [first, *sequence[first]],
                    [last, *sequence[last]],
                    nan_idx
                )
                
                if axis == 0:
                    output_matrix[idx][nan_idx] = rgb # row
                else:
                    output_matrix[:,idx][nan_idx] = rgb # col


def combine_matrices(a, b):
    """Combines two same shaped arrays.
    If there are overlapping RGB values, the average (mean) of both values is used
    """
    assert a.shape == b.shape
    
    t = numpy.array(a, copy=True)
    for x in range(b.shape[0]):
        for y in range(b.shape[1]):
            #print(b[x][y][0], np.isfinite(b[x][y][0]))
            if not numpy.isfinite(b[x][y][0]):
                continue
                
            # If a value at the coords x, y in the origin array is found,
            # use the average of both values
            if numpy.isfinite(a[x][y][0]):
                t[x][y] = numpy.mean((a[x][y], b[x][y]), axis=0) # numpy.average
            # otherwise just set the value of the b array
            else:
                t[x][y] = b[x][y] 
    return t

           
def extrapolate_matrix(matrix):
    
    round_trip = 1
    
    # loop as long as there are empty fields
    while len(numpy.where(numpy.isnan(matrix[:,:,0]))[0]) > 0:
        
        # extrapolate on the y-axis
        exp_y_matrix = numpy.empty(matrix.shape, dtype=matrix.dtype) * nan
        extrapolate_axis(matrix, axis=0, output_matrix=exp_y_matrix)
        if DEBUG:
            add_plot(exp_y_matrix, "Extrapolation on y-axis\n(round: {})".format(round_trip))
        
        # extrapolate on the x-axis
        exp_x_matrix = numpy.empty(matrix.shape, dtype=matrix.dtype) * nan
        extrapolate_axis(matrix, axis=1, output_matrix=exp_x_matrix)
        if DEBUG:
            add_plot(exp_x_matrix, "Extrapolation on x-axis\n(round: {})".format(round_trip))
        
        # combine both matrices of the axes and if necessary calculate the mean of overlapping RGB values
        exp_y_x_matrix = combine_matrices(exp_y_matrix, exp_x_matrix)
        if DEBUG:
            add_plot(exp_y_x_matrix, "Combining extrapolation results\nof x and y-axis (round: {})".format(round_trip))
        
        # combine the original matrix with the exptrapolated one
        matrix = combine_matrices(matrix, exp_y_x_matrix)
        if DEBUG:
            add_plot(matrix, "Combining merged extrapolation\nwith original data (round: {})".format(round_trip))
    
        round_trip += 1
        
    return matrix
    

def resample_matrix(matrix, x_dim, y_dim) :
    return interpolation.zoom(matrix, (x_dim/matrix.shape[0], y_dim/matrix.shape[1], 1), order=1)

   
result_matrix = extrapolate_matrix(filled_grid)

resampled_matrix = resample_matrix(result_matrix, 256, 256)


if DEBUG:
    add_plot(resampled_matrix, "Resample to\n256 colors per axis")
    fig.tight_layout()
    pyplot.show()
else:
    from PIL import Image
    img = Image.fromarray(numpy.array(resampled_matrix, dtype=numpy.uint8), 'RGB')
    img.save('interpolation_result.png')