#!/usr/bin/env python

import re, math, png

path_in          = 'mya/Earthmap1000x500.png'
path_out         = 'Valeriepieris_circle_azimuthal_equal_area.png'
colour_circle    = [255, 255, 0]
radius_circle    = 0.51
thickness_circle = 0.01
lat_centre       = 21.7
long_centre      = 99.383333
zoom             = 0.5
# zoom             = 0.33
# out_size         = 512
out_size         = 2048
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)

## Formula from http://mathworld.wolfram.com/AzimuthalEquidistantProjection.html
def azimuthal_equidistant_to_equirectangular(x, y, lat_centre, long_centre):
 c              = math.hypot(x, y)
 if c == 0 or (abs(lat_centre) == 90 and y == 0): return (0, 0)
 sin_c          = math.sin(c)
 cos_c          = math.cos(c)
 lat_centre_rad = math.radians(lat_centre)
 sin_lat_centre = math.sin(lat_centre_rad)
 cos_lat_centre = math.cos(lat_centre_rad)
 to_asin  = cos_c * sin_lat_centre + y * sin_c * cos_lat_centre / c
 if abs(to_asin) > 1: return (0, 0)
 lat  =  math.degrees(math.asin(to_asin))
 long = (math.degrees(math.atan2(-x, y) if lat_centre ==  90 else
                      math.atan2( x, y) if lat_centre == -90 else
                      math.atan2(x * sin_c, c * cos_lat_centre * cos_c -
                                            y * sin_lat_centre * sin_c)) +
         long_centre + 540) % 360 - 180 ## +540%360-180 to make range [-180,180)
 return (lat, long)
## From http://mathworld.wolfram.com/LambertAzimuthalEqual-AreaProjection.html
def azimuthal_equal_area_to_equirectangular(x, y, lat_centre, long_centre):
 rho            = math.hypot(x, y)
 if rho == 0 or (abs(lat_centre) == 90 and y == 0) or abs(rho * 0.5) > 1:
  return (None, None)
 c              = 2 * math.asin(rho * 0.5)
 sin_c          = math.sin(c)
 cos_c          = math.cos(c)
 lat_centre_rad = math.radians(lat_centre)
 sin_lat_centre = math.sin(lat_centre_rad)
 cos_lat_centre = math.cos(lat_centre_rad)
 to_asin  = cos_c * sin_lat_centre + y * sin_c * cos_lat_centre / rho
 if abs(to_asin) > 1: return (None, None)
 lat  =  math.degrees(math.asin(to_asin))
 long = (math.degrees(math.atan2(x * sin_c, rho * cos_lat_centre * cos_c -
                                            y * sin_lat_centre * sin_c)) +
         long_centre + 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)
   (lat, long) = azimuthal_equal_area_to_equirectangular(x, y, lat_centre, long_centre)
   if lat is None or long is None:
    colour = [0] * png_out.planes
   else:
    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
png_out.write(path_out)