#!/usr/bin/env python

import re, math, png

path_in          = 'Earthmap1000x500.png'
path_out         = 'Mecca_azimuthal_equidistant.png'
colour_circle    = [255, 255, 0]
radius_circle    = 90.51
thickness_circle = 0.01
# lat_centre       = 52.205278
# long_centre      = 0.119167
# lat_centre       = 21.4225
# long_centre      = -1.7582
lat_centre       = 21.3891
long_centre      = 39.8579
zoom             = 0.43
out_size         = 1024
out_size_half    = out_size * 0.5

class Png:
 def __init__(self, path_in):
  (self.width, self.height, self.pixels, self.metadata) = png.Reader(path_in).read_flat()
  self.planes = self.metadata['planes']
 def __str__(self): return str((self.width, self.height, len(self.pixels), self.metadata))
 def write(self, path_out):
  png.Writer(width=self.width, height=self.height,
             bitdepth=self.metadata['bitdepth'], interlace=self.metadata['interlace'],
             planes=self.metadata['planes'], greyscale=self.metadata['greyscale'],
             alpha=self.metadata['alpha']).write_array(open(path_out, 'wb'), self.pixels)

## Formulas from http://mathworld.wolfram.com/AzimuthalEquidistantProjection.html
def azimuthal_equidistant_to_equirectangular(x, y, lat1, long1):
 c        = math.hypot(x, y)
 if c == 0 or (abs(lat1) == 90 and y == 0): return (0, 0)
 sin_c    = math.sin(c)
 cos_c    = math.cos(c)
 lat1_rad = math.radians(lat1)
 sin_lat1 = math.sin(lat1_rad)
 cos_lat1 = math.cos(lat1_rad)
 to_asin  = cos_c * sin_lat1 + y * sin_c * cos_lat1 / c
 if abs(to_asin) > 1: return (0, 0)
 lat  =  math.degrees(math.asin(to_asin))
 long = (math.degrees(math.atan2(-x, y) if lat1 ==  90 else
                      math.atan2( x, y) if lat1 == -90 else
                      math.atan2(x * sin_c, c * cos_lat1 * cos_c - y * sin_lat1 * sin_c)) +
         long1 + 540) % 360 - 180 ## + 540 % 360 - 180 to make range [-180, 180)
 return (lat, long)

png_in  = Png(path_in)
print(png_in)
print(png_in.pixels[:20])
png_out = Png(path_in) ## copy most of original's metadata
png_out.width  = png_out.height = out_size
png_out.pixels = [0] * (png_out.width * png_out.height)
print(png_out)
for  out_y in range(out_size):
 for out_x in range(out_size):
  x = (out_x / out_size_half - 1) /  zoom
  y = (out_y / out_size_half - 1) / -zoom
  if abs(math.hypot(x,y) - radius_circle) < thickness_circle * zoom:
   colour = colour_circle
  else:
   (lat, long) = azimuthal_equidistant_to_equirectangular(x, y, lat_centre, long_centre)
   in_y = int(png_in.height * ( 90 - lat ) / 180.0)
   in_x = int(png_in.width  * (180 + long) / 360.0)
   in_offset = (in_y  * png_in.width + in_x ) * png_in .planes
   colour = png_in.pixels[in_offset :in_offset  + png_in.planes]
  out_offset = (out_y * out_size     + out_x) * png_out.planes
  png_out.pixels[out_offset:out_offset + png_out.planes] = colour
  # if x == 0: print(out_y, out_x, colour)
png_out.write(path_out)