Bitmap_Rectangular_to_BorealTriaxial_4.py

from PIL import Image
from math import atan2, pi, sqrt
import os
import platform

project_folder = 'Maps/TriaxialProjection/'

source_map_path = os.path.join(project_folder, 'Land_shallow_topo_2048.jpg')

output_map_path = os.path.join(project_folder, 'Generated_map_4_4400_shallow.jpg')
output_map_size = 4400
background = (255, 255, 255)
deg = pi/180
stretch = 0.8

sectors = [
    {'interaxial': -22.5 * deg, 'axial': 20 * deg},  # 0
    {'axial': 20 * deg, 'interaxial': 63 * deg},  # 1
    {'interaxial': 63 * deg, 'axial': 130 * deg},  # 2
    {'axial': 130 * deg, 'interaxial': (-150 + 360) * deg},  # 3
    {'interaxial': (-150 + 360) * deg, 'axial': (-70 + 360) * deg},  # 4
    {'axial': (-70 + 360) * deg, 'interaxial': (-22.5 + 360) * deg},  # 5
    {'interaxial': (-22.5 + 360) * deg, 'axial': (20 + 360) * deg},  # 6
]

rotate = 20*deg


def normalize(xy, size):
    xy = list(xy)
    if xy[0] >= size[0]:
        xy[0] = xy[0] - size[0]
    if xy[1] >= size[1]:
        xy[1] = size[1]-1
    return tuple(xy)


def project_point(i, j, source_map, output_map, output_size):
    x = (i + 0.5) - output_size / 2
    y = (j + 0.5) - output_size / 2
    a_p = (atan2(x, y) + rotate) % (2 * pi)
    r = sqrt(x * x + y * y) / (output_size / 2) * pi
    for s in sectors:
        if s['range'][0] <= a_p <= s['range'][1]:
            axial = s['axial']
            interaxial = s['interaxial']
    f_p = (a_p - interaxial) / (axial - interaxial)
    if r < pi:
        f_r = f_p ** (stretch*pi / (pi - r))
    else:
        f_r = 0
    a_r = (axial - interaxial) * f_r + interaxial
    lat = pi / 2 - r
    lon = (a_r + pi) % (2*pi) - pi
    point = (
        (lon + pi) / (2 * pi) * source_map.size[0],
        (pi / 2 - lat) / pi * source_map.size[1],
    )
    neighbors = [
        {'coord': (int(point[0]) + 0, int(point[1]) + 0), 'close': (1 - point[0] % 1) * (1 - point[1] % 1)},
        {'coord': (int(point[0]) + 0, int(point[1]) + 1), 'close': (1 - point[0] % 1) * (0 + point[1] % 1)},
        {'coord': (int(point[0]) + 1, int(point[1]) + 0), 'close': (0 + point[0] % 1) * (1 - point[1] % 1)},
        {'coord': (int(point[0]) + 1, int(point[1]) + 1), 'close': (0 + point[0] % 1) * (0 + point[1] % 1)},
    ]
    for n in neighbors:
        n['coord'] = normalize(n['coord'], source_map.size)
        n['color'] = source_map.getpixel(n['coord'])
    color = tuple([int(sum([n['color'][channel]*n['close'] for n in neighbors])) for channel in range(3)])

    if point[0] < source_map.size[0] and point[1] < source_map.size[1]:
        output_map.putpixel((i, j), color)


def project_rectangular_bitmap_to_boreal_triaxial(source=source_map_path, output_size=360):
    for s in sectors:
        s['range'] = sorted([s['axial'], s['interaxial']])
    source_map = Image.open(source).convert('RGB')
    r, g, b = source_map.getpixel((0, 0))
    output_map = Image.new('RGB', (output_size, output_size), background)
    for i in range(output_size):
        for j in range(output_size):
            project_point(i, j, source_map, output_map, output_size)
    output_map.save(output_map_path)


project_rectangular_bitmap_to_boreal_triaxial(source_map_path, output_map_size)