A Python script I hacked together to convert OS mapping data to Minecraft maps.

osmap-to-mc.py

#!/usr/bin/env python

from __future__ import print_function

import os
import numpy as np
import glob

import mclevel
import materials

import gdal

terrain_data_root = '/media/rjw57/Celsus/os-terrain-50/data'
colour_data_root = '/media/rjw57/Celsus/os-miniscale'
level = mclevel.fromFile(os.path.expanduser('~/.minecraft/saves/os-high-real-height/'))

terrain_data_files = glob.glob(os.path.join(terrain_data_root, '*', '*.asc'))
colour_data_files = glob.glob(os.path.join(colour_data_root, 'MiniScale_relief1_R15.tif'))

SEA_HEIGHT = 63
BLOCK_SIZE = 100
BLOCK_HEIGHT = 25

def get_or_create_chunk(level, cx, cy):
    try:
        return level.getChunk(cx, cy)
    except:
        level.createChunk(cx, cy)
        c = level.getChunk(cx, cy)
        c.Blocks[:,:,:] = materials.indevMaterials.Air.ID
        c.Blocks[:,:,:SEA_HEIGHT] = materials.indevMaterials.Dirt.ID
        c.Blocks[:,:,SEA_HEIGHT] = materials.indevMaterials.Water.ID
        c.chunkChanged()
        return c

def process_elevation_chunk(level, cx, cy, block_to_pixel, data, height_scale):
    c = get_or_create_chunk(level, cx, cy)
    for bx in xrange(16):
        x = (cx<<4) + bx
        for by in xrange(16):
            y = (cy<<4) + by
            
            px, py, _ = block_to_pixel.dot((x, y, 1))
            
            if px < 0 or px >= data.shape[1]:
                continue
            if py < 0 or py >= data.shape[0]:
                continue
            
            h_m = data[py, px]
            h = np.clip(np.floor(SEA_HEIGHT + h_m*height_scale), 0, 255)
            c.Blocks[bx,by,h:] = materials.indevMaterials.Air.ID
            
            if h_m < 0:
                c.Blocks[bx,by,:SEA_HEIGHT+1] = materials.indevMaterials.Water.ID
                c.Blocks[bx,by,:h+1] = materials.indevMaterials.Dirt.ID
            elif h_m < 500:
                c.Blocks[bx,by,:h+1] = materials.indevMaterials.Grass.ID
            elif h_m < 1000:
                c.Blocks[bx,by,:h+1] = materials.indevMaterials.Stone.ID
            else:
                c.Blocks[bx,by,:h+1] = 80 # Snow

            c.Blocks[:,:,0] = materials.indevMaterials.Bedrock.ID
            
    c.Biomes[:,:] = 1 # Plains
    c.chunkChanged()

def load_tile(tile_filename, block_size):
    tile = gdal.Open(tile_filename)

    # Create a matrix which maps pixel coords -> projection
    ox, sx, _, oy, _, sy = tile.GetGeoTransform()
    pix_to_proj = np.array(((sx, 0, ox), (0, sy, oy), (0,0,1)))
    
    proj_to_pix = np.linalg.inv(pix_to_proj)
    
    # How big is one minecraft block in projection coords?
    block_to_proj = np.diag((block_size, -block_size, 1))
    
    # Create a matrix which maps minecraft block coords -> pixel coords
    # this is block -> proj -> pixel
    block_to_pixel = proj_to_pix.dot(block_to_proj)
    pixel_to_block = np.linalg.inv(block_to_pixel)
    
    # Compute block extent
    x1, y1, _ = pixel_to_block.dot((0,0,1))
    x2, y2, _ = pixel_to_block.dot((tile.RasterXSize,tile.RasterYSize,1))
    minx = int(min(x1, x2))
    miny = int(min(y1, y2))
    maxx = int(max(x1, x2))
    maxy = int(max(y1, y2))
    
    # Read tile data
    array = tile.ReadAsArray()

    return array, block_to_pixel, (minx, maxx, miny, maxy)

def process_elevation_tile(level, tile_filename):
    dtm_array, block_to_pixel, extent = load_tile(tile_filename, block_size=BLOCK_SIZE)
    
    minx, maxx, miny, maxy = extent
    for x in xrange(minx>>4, ((maxx-1)>>4)+1):
        for y in xrange(miny>>4, ((maxy-1)>>4)+1):
            process_elevation_chunk(level, x, y, block_to_pixel, dtm_array, 1.0/BLOCK_HEIGHT)
            
    return extent

def process_colour_chunk(level, cx, cy, block_to_pixel, data):
    c = get_or_create_chunk(level, cx, cy)
    for bx in xrange(16):
        x = (cx<<4) + bx
        for by in xrange(16):
            y = (cy<<4) + by
            
            px, py, _ = block_to_pixel.dot((x, y, 1))
            
            if px < 0 or px >= data.shape[1]:
                continue
            if py < 0 or py >= data.shape[0]:
                continue
            
            colour = data[py, px, :]
            wool_id = best_wool(colour)

            column = c.Blocks[bx,by,:]
            to_colour = np.logical_and(column != materials.indevMaterials.Air.ID, column != materials.indevMaterials.Water.ID)

            c.Blocks[bx,by,to_colour] = 35 # wool
            c.Data[bx,by,to_colour] = wool_id
            
    c.Biomes[:,:] = 1 # Plains
    c.chunkChanged()

def process_colour_tile(level, tile_filename):
    colour_array, block_to_pixel, extent = load_tile(tile_filename, block_size=BLOCK_SIZE)

    # The colour array has shape 3xheightxwidth which is not quite convenient
    colour_array = np.transpose(colour_array, (1,2,0))
    
    minx, maxx, miny, maxy = extent
    last_progess = -np.inf
    for x in xrange(minx>>4, ((maxx-1)>>4)+1):
        progress = 100.0 * (x-(minx>>4)) / float(maxx>>4)
        if progress >= last_progess + 1:
            print('Colour tile progress: {0}%'.format(int(progress)))
            last_progess = progress

        for y in xrange(miny>>4, ((maxy-1)>>4)+1):
            process_colour_chunk(level, x, y, block_to_pixel, colour_array)
            
    return extent

def hex_to_rgb(hex_code):
    """Convert a hex value of the form 0xXXRRGGBB to a RGB colour value."""
    return np.array(((hex_code >> 16) & 0xff, (hex_code >> 8) & 0xff, hex_code & 0xff))

# These are the approximate wool colours ordered by datavalue
WOOL_COLOURS = np.array([
    hex_to_rgb(0xe4e4e4),
    hex_to_rgb(0xea7e35),
    hex_to_rgb(0xbe49c9),
    hex_to_rgb(0x6387d2),
    hex_to_rgb(0xc2b51c),
    hex_to_rgb(0x39ba2e),
    hex_to_rgb(0xd98199),
    hex_to_rgb(0x414141),
    hex_to_rgb(0xa0a7a7),
    hex_to_rgb(0x267191),
    hex_to_rgb(0x7e34bf),
    hex_to_rgb(0x253193),
    hex_to_rgb(0x56331c),
    hex_to_rgb(0x364b18),
    hex_to_rgb(0x9e2b27),
    hex_to_rgb(0x181414),
], dtype=np.float32)

def best_wool(colour):
    """Given a RGB colour triplet, find the best matching colour."""

    # Compute a 3xN array of differences
    deltas = np.array(list(WOOL_COLOURS[:,i] - colour[i] for i in xrange(3)))

    # Compute squared distance
    sq_distances = np.sum(deltas*deltas, axis=0)

    return np.argmin(sq_distances)

print('Processing terrain')
files = terrain_data_files
last_progess = -np.inf
extent = (np.inf, -np.inf, np.inf, -np.inf)
for idx, fn in enumerate(files):
    progress = 100.0 * (idx+1.0) / len(files)
    if progress >= last_progess+1:
        print('{0}%'.format(int(progress)))
        last_progess = progress
        
    minx, maxx, miny, maxy = process_elevation_tile(level, fn)
    extent = (
        min(extent[0], minx), max(extent[1], maxx),
        min(extent[2], miny), max(extent[3], maxy)
    )
print('Terrain Extent: {0}'.format(extent))   

print('Processing colour')
files = colour_data_files
extent = (np.inf, -np.inf, np.inf, -np.inf)
last_progess = -np.inf
for idx, fn in enumerate(files):
    progress = 100.0 * (idx+1.0) / len(files)
    if progress >= last_progess+1:
        print('{0}%'.format(int(progress)))
        last_progess = progress
        
    minx, maxx, miny, maxy = process_colour_tile(level, fn)
    extent = (
        min(extent[0], minx), max(extent[1], maxx),
        min(extent[2], miny), max(extent[3], maxy)
    )
print('Colour Extent: {0}'.format(extent))   

level.setPlayerPosition((0.5*(extent[0]+extent[1]), 128, 0.5*(extent[2]+extent[3])))
level.setPlayerSpawnPosition((0.5*(extent[0]+extent[1]), 128, 0.5*(extent[2]+extent[3])))
level.saveInPlace()

# vim:sw=4:sts=4:et