albaenia.py

# coding: utf-8

# In[147]:

from Circles.circles import circle
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
from matplotlib.collections import PatchCollection
from shapely.geometry import Polygon, MultiPolygon, Point, LineString
from descartes import PolygonPatch
import numpy as np
import math
get_ipython().magic(u'matplotlib inline')


# In[97]:

# won't work unless you have LaTeX and Helvetica
from matplotlib import rc
rc('font', **{'family':'sans-serif',
    'sans-serif':['Helvetica'],
    'monospace': ['Inconsolata'],
    'serif': ['Helvetica']})
rc('text', **{'usetex': True})
rc('text', **{'latex.preamble': '\usepackage{sfmath}'})


# In[196]:

centerlon = -0.190278
centerlat = 51.148056


# In[192]:

def triangle_area(*vertices):
    """
    Use Heron's formula to calculate area of a triangle
    Returns answer in km^2, assuming input in m
    """
    s = (vertices[0].length + vertices[1].length + vertices[2].length) / 2
    return math.sqrt(s * (s - vertices[0].length) * (s - vertices[1].length) * (s - vertices[2].length)) / 1000000


# Great-circle lines drawn on a Lambert Conformal projection approximate straight lines for flight distances

# In[217]:

# setup lambert conformal basemap.
# lat_1 is first standard parallel.
# lat_2 is second standard parallel (defaults to lat_1).
# lon_0, lat_0 is central point.
# rsphere=(6378137.00,6356752.3142) specifies WGS4 ellipsoid
# area_thresh=1000 means don't plot coastline features less
# than 1000 km^2 in area.
m = Basemap(width=6000000, height=5000000,
            rsphere=(6378137.00, 6356752.3142),
            resolution='h', area_thresh=1000., projection='lcc',
            lat_1=51., lat_2=51,
            lat_0=centerlat,lon_0=centerlon)


# In[218]:

RADIUS = 50
coords = circle(m, centerlon, centerlat, RADIUS)
home = Polygon(coords)


# In[219]:

# Tirana
tlon = 19.831799999999930000
tlat = 41.331650000000000000
tirana = Polygon(circle(m, tlon, tlat, RADIUS))


# In[220]:

# Tunis
tunlon = 10.165960000000041000
tunlat = 36.818810000000000000
tunis = Polygon(circle(m, tunlon, tunlat, RADIUS))


# In[228]:

# triangle
area = Polygon((m(centerlon, centerlat), m(tlon, tlat), m(tunlon, tunlat)))


# In[242]:

# let's plot everything
plt.clf()
fig = plt.figure()
ax = fig.add_subplot(111, axisbg='w', frame_on=False)

home_patch = PolygonPatch(home, fc='#CC00CC', ec='#333333', lw=1.5, alpha=.75, zorder=2)
tir_patch = PolygonPatch(tirana, fc='#CC00CC', ec='#333333', lw=1.5, alpha=.75, zorder=2)
tun_patch = PolygonPatch(tunis, fc='#CC00CC', ec='#333333', lw=1.5, alpha=.75, zorder=2)
area_patch = PolygonPatch(area, fc='none', ls='dotted', hatch='...', ec='none', alpha=.75, zorder=2)

# draw map features
m.drawmapboundary(fill_color='#dfedff', linewidth=0.25, zorder=0)
# m.drawcountries(linewidth=.25, color='#000000', zorder=2)
m.drawcoastlines(linewidth=1., zorder=1)
m.fillcontinents(color='#555555', lake_color='#555555',zorder=1)
m.drawparallels(np.arange(-90., 120., 30.), alpha=0.5, lw=0.25, zorder=1)
m.drawmeridians(np.arange(0., 360., 60.), alpha=0.5, lw=0.25, zorder=1)

# draw circles + triangle
ax.add_patch(home_patch)
ax.add_patch(tir_patch)
ax.add_patch(tun_patch)
ax.add_patch(area_patch)

lon_gc = m.drawgreatcircle(centerlon, centerlat, tlon, tlat, lw=1.75, color='#FFCC00', zorder=1)
tir_gc = m.drawgreatcircle(tlon, tlat, tunlon, tunlat, lw=1.75, color='#FFCC00', zorder=1)
tun_gc = m.drawgreatcircle(tunlon, tunlat, centerlon, centerlat, lw=1.75, color='#FFCC00', zorder=1)

distances = [LineString(leg) for leg in [lon_gc._vertices, tir_gc._vertices, tun_gc._vertices]]

m.drawmapscale(
    m.llcrnrlon + 1.75, m.llcrnrlat + 3.25,
    centerlon, centerlat,
    500.,
    barstyle='fancy', labelstyle='simple',
    fillcolor1='w', fillcolor2='#555555',
    fontcolor='#555555', fontsize=12,
    zorder=5)

plt.title(
    r"LGW $\rightarrow$ TIA, TIA $\rightarrow$ TUN, TUN $\rightarrow$ LGW. %s$km^2$ of Misery" % "{:.2f}".format(triangle_area(*distances)),
    size=16)
plt.tight_layout()
fig.set_size_inches(12., 12.)
plt.savefig('albaenia.png', dpi=300, bbox_inches='tight', alpha=False, transparent=False)
plt.show()


# In[ ]: