I wanted to plot Mexico and Switzerland on top of each other, with the same scale, to compare the size and relative distances. (If you try that in Google maps you run into troubles, as the scale is not the same; try it!)
There is nowadays a wealth of data in the Open Domain when it comes to administrative borders. For this example I use the following data:
More resources can be found on the wiki of the Open Street Map, the U.S. Humanitarian Information Unit [HIU], and an internet search will provide many more.
In [1]:
import shapefile
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
Settings for both figures
I want that Mexico City and my home village Thalheim overlay (I took the coordinates from Wikipedia). The easiest way to achieve this is to use them as the centre of the map. However, I shift them subsequently by 2° to the south, as Mexico City is not in the centre of Mexico.
In [2]:
# Coordinates (from Wikipedia) and shift
DF = [-99.1333, 19.4333]
TH = [ 8.1042, 47.4375]
shift = [0, 2]
# Colours
cMX = '#0000cc' # Blue tone
cCH = '#006600' # Green tone
Figure 1.A: Transverse Mercator
In [3]:
# 4000 km east-west, 3000 km north-south, low resolution, Transverse Mercator
width = 4000000
height = 3000000
res = 'l'
proj = 'tmerc'
fig1a = plt.figure(figsize=(8, 6))
# Create basemaps for Mexico and Switzerland
m_mx = Basemap(width=width, height=height, resolution=res, projection=proj,
lon_0=DF[0]+shift[0], lat_0=DF[1]+shift[1])
m_ch = Basemap(width=width, height=height, resolution=res, projection=proj,
lon_0=TH[0]+shift[0], lat_0=TH[1]+shift[1])
# Draw coast lines
m_mx.drawcoastlines(color=cMX, linewidth=1)
m_ch.drawcoastlines(color=cCH, linewidth=1)
# Draw all country borders
m_mx.drawcountries(color=cMX, linewidth=1)
m_ch.drawcountries(color=cCH, linewidth=1)
# Fill continents and lakes
m_mx.fillcontinents(lake_color='none', color=cMX, alpha=.4)
m_ch.fillcontinents(lake_color='none', color=cCH, alpha=.4)
# Draw extra points for the Swiss and Mexican countries to highlight
CHE_adm0 = shapefile.Reader('data/basemap/CHE_adm0')
MEX_adm0 = shapefile.Reader('data/basemap/MEX_adm0')
chlonlat = np.array(CHE_adm0.shape().points)
mxlonlat = np.array(MEX_adm0.shape().points)
m_ch.plot(chlonlat[:, 0], chlonlat[:, 1] , '.', c=cCH, latlon=True, ms=1)
m_mx.plot(mxlonlat[:, 0], mxlonlat[:, 1], '.', c=cMX, latlon=True, ms=1)
# Draw scales to cross-check that the scales are the same
# (shift scale 1600 km to the west from centre, and a bit north/south)
smx = np.array(m_mx(DF[0], DF[1])) - np.array([1600000, 50000])
sch = np.array(m_ch(TH[0], TH[1])) - np.array([1600000, -350000])
ismx = m_mx(smx[0], smx[1], inverse='True')
isch = m_ch(sch[0], sch[1], inverse='True')
m_mx.drawmapscale(ismx[0], ismx[1], DF[0]+shift[0], DF[1]+shift[1],
500, barstyle='fancy', fillcolor2=cMX, fontsize=12)
m_ch.drawmapscale(isch[0], isch[1], TH[0]+shift[0], TH[1]+shift[1],
500, barstyle='fancy', fillcolor2=cCH, fontsize=12)
# Add Mexico City and Thalheim AG
m_mx.plot(DF[0], DF[1], 'x', c='r', ms=10, mew=2, latlon=True, label='Thalheim')
m_ch.plot(TH[0], TH[1], '+', c='r', ms=14, mew=2, latlon=True, label='Mexico City')
lg = plt.legend(loc='lower left', fontsize=12, numpoints=1)
lg.get_frame().set_alpha(.8) # A little transparency
# Draw a thick border around the whole
m_mx.drawmapboundary(color='k', linewidth=3)
plt.show()
Figure 1.B: Orthographic
In [4]:
fig1b = plt.figure(figsize=(8, 8))
# Create basemaps for Mexico and Switzerland
m_mx = Basemap(resolution='c', projection='ortho', lon_0=DF[0], lat_0=DF[1])
m_ch = Basemap(resolution='c', projection='ortho', lon_0=TH[0], lat_0=TH[1])
# Draw coast lines
m_mx.drawcoastlines(color=cMX, linewidth=1)
m_ch.drawcoastlines(color=cCH, linewidth=1)
# Draw all country borders
m_mx.drawcountries(color=cMX, linewidth=1)
m_ch.drawcountries(color=cCH, linewidth=1)
# Fill continents and lakes
m_mx.fillcontinents(lake_color='none', color=cMX, alpha=.4)
m_ch.fillcontinents(lake_color='none', color=cCH, alpha=.4)
plt.show()
Figure 2
In [5]:
# 70 km by 70 km, Transverse Mercator, resolution does not
# matter, as we use other data for the borders
width = 70000
height = 70000
dshift = [-.133, -.05]
fig2 = plt.figure(figsize=(8,8))
# Create basemaps for Mexico and Switzerland
m_mx = Basemap(width=width, height=height, projection=proj,
lon_0=DF[0], lat_0=DF[1]+dshift[0])
m_ch = Basemap(width=width, height=height, projection=proj,
lon_0=TH[0], lat_0=TH[1]+dshift[1])
# Draw the district of Aargau
CHE_adm1 = shapefile.Reader('data/basemap/CHE_adm1')
iag = [i for i, s in enumerate(CHE_adm1.records()) if 'Aargau' in s][0]
aglonlat = np.array(CHE_adm1.shapes()[iag].points)
m_ch.plot(aglonlat[:, 0], aglonlat[:, 1], '-', c=cCH, latlon=True, lw=2)
agx, agy = m_ch(aglonlat[:, 0], aglonlat[:, 1])
plt.fill(agx, agy, cCH, ec='none', alpha=.4)
# Draw Mexico City DF
MEX_adm1 = shapefile.Reader('data/basemap/ne_10m_admin_1_states_provinces')
idf = [i for i, s in enumerate(MEX_adm1.records()) if 'MX.DF' in s][0]
dflonlat = np.array(MEX_adm1.shape(idf).points)
m_mx.plot(dflonlat[:, 0], dflonlat[:, 1], '-', c=cMX, latlon=True, lw=2)
dfx, dfy = m_mx(dflonlat[:, 0], dflonlat[:, 1])
plt.fill(dfx, dfy, cMX, ec='none', alpha=.4)
# Draw scales to cross-check that the scales are the same
sdf = np.array(m_mx(DF[0], DF[1])) - np.array([28000, 10000])
sth = np.array(m_ch(TH[0], TH[1])) - np.array([28000, 10000])
isdf = m_mx(sdf[0], sdf[1], inverse='True')
isth = m_ch(sth[0], sth[1], inverse='True')
m_mx.drawmapscale(isdf[0], isdf[1], DF[0], DF[1]+dshift[0],
4, barstyle='fancy', fillcolor2=cMX, fontsize=12)
m_ch.drawmapscale(isth[0], isth[1], TH[0], TH[1]+dshift[1],
4, barstyle='fancy', fillcolor2=cCH, fontsize=12)
# Draw locations
hPT = [-99.12, 19.47]
pPT = [-99.11, 19.33]
m_mx.plot(hPT[0], hPT[1], '*', c='k', ms=14, mew=1, latlon=True)
m_mx.plot(pPT[0], pPT[1], 'o', c='k', ms=10, mew=1, latlon=True)
# Remove border around the whole
m_mx.drawmapboundary(color='none')
plt.show()
Mexico City (blue) and my home district Aargau (green) are approximately the same size (Wikipedia: DF ~1485 km2, AG ~1404 km2). True, roughly a third to the south and the south-west of the above outline of Mexico City is mountainous area. But then the city extends to the east and the north outside the official DF-boundaries, so the area-comparison roughly holds.
In Mexico, when we go from home (*) to my parents in law (o), we make about the same distance as if we would go from Thalheim to the German border, just without ever leaving the city. Interesting indeed.