ozi2geotiff.py

################################################################################
## ozi2geotiff ( convert maps from Ozi Explorer format to GeoTiff )
##
## Copyright (C) 2008 Andrew Vagin
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2 of the License, or
## (at your option) any later version.
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program; if not, write to the Free Software
## Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA
##
## Andrew Vagin <avagin@gmail.com>
##
################################################################################

from logging import critical, error, warn, info, debug
import os, logging, ntpath
from optparse import OptionParser
from subprocess import Popen, PIPE
import math
import tempfile
DEG_TO_RAD = math.pi / 180.0

def degmin_to_deg(s, d, m):
	deg = abs(d) + m / 60.0;
	if s:
		deg = -deg
	return deg

def gdalwarp(infile, outfile, order = 0):
	info("execute gdalwarp")
	#TODO investigate -order and -tps
	args = ["gdalwarp", "-tps"]
	args.extend(["-dstnodata", "255"])
	args.extend([infile, outfile])
	debug(args)
	gdalwarp = Popen(args)
	ret = gdalwarp.wait()
	if ret:
		raise Exception("gdalwarp failed: return %d" % ret)

def gdal_translate2(infile, outfile):
	info("execute gdal_translate")
	args = ["gdal_translate"]
	args += ["-co", "tiled=yes"]
	args += ["-co", "blockxsize=256"]
	args += ["-co", "blockysize=256"]
	args += ["-co", "compress=deflate"]
	args += ["-co", "predictor=1"]
	args += [infile, outfile]
	debug(args)
	gdal_translate = Popen(args)
	ret = gdal_translate.wait()
	if ret:
		raise Exception("gdal_translate failed: return %d" % ret)

def convert_ozi_map(mapfile, imgfile, outfile):
	# defule projection.
	# TODO: get projection from ozi files
	projection = "+proj=merc +ellps=WGS84 +datum=WGS84 +no_defs"
	fd = open(mapfile)
	points = []
	proj_fd = Popen("proj %s" % projection, \
			stdin = PIPE, stdout = PIPE, shell = True)

	points = []
	i = 0;
	for line in fd:
		i += 1
		if not imgfile and i == 2:
			imgfile = ntpath.basename(line.strip())
			imgfile = os.path.join(os.path.dirname(mapfile), imgfile)
			info("image file: %s" % imgfile)
			if not os.access(imgfile, os.F_OK):
				raise Exception("can't find img file")
		if not line.startswith("Point"):
			continue
		d = line.split(",")
		for i in range(len(d)):
			d[i] = d[i].strip()
		# map file may contain stupid lines
		# Point10
		# Point10,xy,,,,,,
		# and etc
		if len(d) < 3 or not d[2]:
			continue
		x = d[2]
		y = d[3]
		points.append((x,y))
		signLat = d[8]
		degLat = d[6]
		minLat = d[7]
		signLon = d[11]
		degLon = d[9]
		minLon = d[10]
		lat = "%sd%s%s" % (degLat, minLat, signLat)
		lon = "%sd%s%s" % (degLon, minLon, signLon)
		s = "%s %s\n" % (lon, lat)
		debug(s.strip())
		proj_fd.stdin.write(s)
	proj_fd.stdin.close()

	info("convert points")
	args = ["-a_srs", projection];
	for line in proj_fd.stdout:
		debug(line.strip())
		lon, lat = line.split()
		p = points.pop(0)
		args.extend(["-gcp", p[0], p[1], lon, lat])
	ret = proj_fd.wait()
	if ret:
		raise Exception("proj failed: return %d" % ret)

	tmp1 = tempfile.mktemp(suffix=".tif", prefix = "mapconverter_")
	try:
		info("execute gdal_translate")
		args.extend([imgfile, tmp1])
		args = ["gdal_translate"] + args
		debug(args)
		gdal_translate = Popen(args)
		ret = gdal_translate.wait()
		if ret:
			raise Exception("gdal_translate failed: return %d" % ret)
		
		tmp2 = tempfile.mktemp(suffix=".tif", prefix = "mapconverter_")
		tmp3 = tempfile.mktemp(suffix=".tif", prefix = "mapconverter_")
		tmp4 = tempfile.mktemp(suffix=".tif", prefix = "mapconverter_")
		try:
			gdalwarp(tmp1, tmp2)

			gdal_translate2(tmp2, tmp3)
		finally:
			if os.access(tmp2, os.F_OK):
				debug("delete file %s" % tmp2)
				os.unlink(tmp2)
		cmd = ["convert", "-depth", "8", "-type", "Palette", tmp3, tmp4]
		debug(cmd)
		fd = Popen(cmd)
		ret = fd.wait()
		if ret:
			raise Exception("convert failed. exit code = %d" % ret)
		tmp5 = tempfile.mktemp(suffix=".data", prefix = "mapconverter_")
		f = open(tmp5, "w")
		cmd = ["listgeo", "-no_norm", tmp3]
		debug(cmd)
		fd = Popen(cmd, stdout = f)
		ret = fd.wait()
		if ret:
			raise Exception("listgeo failed. exit code = %d" % ret)
		f.close()
		cmd = ["geotifcp", "-g", tmp5, tmp4, outfile]
		debug(cmd)
		fd = Popen(cmd)
		ret = fd.wait()
		if ret:
			raise Exception("geotifcp failed. exit code = %d" % ret)


	finally:
		if os.access(tmp1, os.F_OK):
			debug("delete file %s" % tmp1)
			os.unlink(tmp1)
		if os.access(tmp3, os.F_OK):
			debug("delete file %s" % tmp3)
			os.unlink(tmp3)
		if os.access(tmp4, os.F_OK):
			debug("delete file %s" % tmp4)
			os.unlink(tmp4)
		if os.access(tmp5, os.F_OK):
			debug("delete file %s" % tmp5)
			os.unlink(tmp5)

parser = OptionParser()

parser.add_option('-v', '--verbose', action="store_true",
		help = 'print debug info')
parser.add_option('-o', '--out-file', default = None, \
		help = "result geotiff map",
		type = 'string')
parser.add_option('-i', '--in-file', default = None, \
		help = "converted map file.",
		type = 'string')
parser.add_option('--in-img', default = None, \
		help = "map image file",
		type = 'string')


options, args = parser.parse_args()
if not options.in_file:
	parser.print_help()
	parser.error( ("The options -i is obligatory!") )

if options.verbose:
	logging.getLogger().setLevel(logging.DEBUG)
else:
	logging.getLogger().setLevel(logging.INFO)

convert_ozi_map(options.in_file, options.in_img, options.out_file)