Quick and dirty implementation of the sinusoidal projection.

main.py

#!/usr/bin/env python


import sys
import math
import itertools as it
import numpy as np
import scipy.misc


def latlon_to_xy(lat, lon, central_meridian=0):
	x = (lon - central_meridian) * math.cos(lat)
	y = lat
	return x, y


def xy_to_latlon(x, y, central_meridian=0):
	lat = y
	lon = central_meridian + x / math.cos(lat)
	return lat, lon


def ij_to_latlon(i, j, w, h):
	lon = 2 * math.pi * (j - h/2) / h
	lat = math.pi * (i - w/2) / w
	return lat, lon


def latlon_to_ij(lat, lon, w, h):
	i = lat * w / math.pi + w/2
	j = lon * h / (2 * math.pi) + h/2
	return i, j


def xy_to_ij(x, y, x_min, x_max, y_min, y_max, w, h):
	j = h * (x + x_max) / (x_max - x_min)
	i = w * (y + y_max) / (y_max - y_min)
	return i, j


def ij_to_xy(i, j, x_min, x_max, y_min, y_max, w, h):
	x = (x_max - x_min) * j / h - x_max
	y = (y_max - y_min) * i / w - y_max
	return x, y


def copy_fwd(src, dst, x_min, x_max, y_min, y_max):
	w, h, _ = src.shape
	for i, j in it.product(range(w), range(h)):
		lat, lon = ij_to_latlon(i, j, w, h)
		x, y = latlon_to_xy(lat, lon) 
		I, J = xy_to_ij(x, y, x_min, x_max, y_min, y_max, w, h)
		dst[int(I), int(J), :] = src[i, j, :] 


def copy_bwd(src, dst, x_min, x_max, y_min, y_max):
	w, h, _ = src.shape
	for I, J in it.product(range(w), range(h)):
		x, y = ij_to_xy(I, J, x_min, x_max, y_min, y_max, w, h)
		lat, lon = xy_to_latlon(x, y)
		i, j = latlon_to_ij(lat, lon, w, h)
		if min(i, j) >= 0 and i < w and j < h:
			dst[I, J, :] = src[int(i), int(j), :]


if __name__ == '__main__':
	if len(sys.argv) <= 2:
		print('Usage: python ', sys.argv[0], ' input_image output_image', sep=' ')
		sys.exit(1)

	in_fname = sys.argv[1]
	out_fname = sys.argv[2]

	im_in = scipy.misc.imread(in_fname, mode="RGBA")

	w, h, d = im_in.shape
	print('w=', w, 'h=', h, 'd=', d)
	im_out = np.zeros([w, h, d])
	im_out[:, :, 3] = 1

	x_max, y_max = math.pi, math.pi / 2
	x_min, y_min = -math.pi, -math.pi / 2

	print('x_min=', x_min, 'x_max=', x_max)
	print('y_min=', y_min, 'y_max=', y_max)
	print(x_max - x_min)

	copy_bwd(im_in, im_out, x_min, x_max, y_min, y_max)

	scipy.misc.imsave(out_fname, im_out)