import numpy as np
from pixell import enmap, utils, wcsutils, enplot
from PIL import Image

bg  = np.array(Image.open("background.png"))[:,::-1]
wcs = wcsutils.projection("mer")
shape, wcs = wcsutils.pixelization(wcs, shape=(bg.shape[0],bg.shape[0]))
shape = (shape[0],bg.shape[1])
wcs.wcs.crpix[0] = 485.5
# Find area < 5000 km away from Paris
r    = 5000e3
ang  = r/utils.R_earth
pos  = np.array([48.86, 2.36])*utils.degree
# Something goes wrong with distance_from here
#mask = enmap.distance_from(shape, wcs, pos, rmax=ang)<=ang
dist = enmap.enmap(utils.angdist(enmap.posmap(shape, wcs)[::-1], pos[:,None,None][::-1]),wcs)*utils.R_earth
mask = dist<=r
# Build plot
lmask = enplot.plot(mask, grid=False, mask=0, min=0, max=1, slice="[:,::-1]", color="0:ff000080,1:ff000080", layers=True)[0]
lcont = enplot.plot(dist, grid=False, mask=0, min=0, max=r, contours="1000e3", contour_width=2, slice="[:,::-1]", layers=True)[-1]
lmid  = enplot.plot(dist, grid=False, mask=0, min=0, max=r, contours="0,10e3", contour_width=2, slice="[:,::-1]", layers=True)[-1]
lbg   = enplot.plot(enmap.enmap(np.moveaxis(bg,2,0),wcs), min=0, max=255, rgb=True, grid=False, slice="[:,::-1]",layers=True)[0]
plot  = enplot.merge_plots([lbg,lmask,lcont,lmid])

enplot.write("paris_mercator.png", plot)