Composite all hillshade layers into a single GeoTIFF

merge_mapnik.rb

	# Combine all hillshade layers into a single GeoTIFF
	# Requires gdalcopyproj.py (which uses GDAL Python bindings)

	require '/usr/local/share/ruby/ffi-mapnik/lib/ffi-mapnik.rb'
	require 'nokogiri'

	Mapnik.register_datasources("/usr/local/lib/mapnik/input")
	Mapnik.register_fonts("/usr/local/share/maps/style/fonts")

	# Choose directory and stylesheet paths
	dir = "/usr/local/share/maps/relief"
	ss  = "/usr/local/share/maps/style/stylesheet.xml"

	# Create a copy of the Mapnik style with all non-hillshade layers disabled
	doc = Nokogiri::XML(File.open(ss))
	doc.search("Layer").each do |layer|
		next if ["hillshade","relief","coloured_relief"].include?(layer['name'])
		layer['status'] = 'off'
	end
	File.write("/tmp/hillshade.xml", doc.to_xml)

	# Iterate through all files...
	Dir.glob("#{dir}/*_hs3.tif").each do |fn|
		puts File.basename(fn)
		outfn = fn.gsub("_hs3.tif","_all.tif")

		# Read image dimensions and coordinate information from original GeoTIFF
		gdal = `gdalinfo "#{fn}"`
		if gdal=~/Size is (\d+), (\d+)/ then width=$1.to_i; height=$2.to_i else puts "Couldn't find size"; exit end
		if gdal=~/Upper Left  .\s*([\-\d\.]+),\s*([\-\d\.]+)/ then w=$1.to_f; n=$2.to_f else puts "Couldn't find upper left "; exit end
		if gdal=~/Lower Right .\s*([\-\d\.]+),\s*([\-\d\.]+)/ then e=$1.to_f; s=$2.to_f else puts "Couldn't find lower right"; exit end
		puts "-- size #{width}x#{height}, bounds #{w},#{s},#{e},#{n}"

		# Generate image with Mapnik
		map = Mapnik::Map.new(10,10)
		outbox = Mapnik::Bounds.new(w,s,e,n)
		map.load("/tmp/hillshade.xml")
		map.resize(width,height)
		map.zoom_to(outbox)
		map.to_file(outfn)
		map.free
		outbox.free

		# Reinstate original GeoTIFF information
		`gdalcopyproj.py "#{fn}" "#{outfn}"`
		puts "-> #{outfn}"
	end