With global warming melting the icecaps and opening up the poles for oil exploration and tourism, I think it’s time for a new standard wall map, one that shifts those distorted map regions away from major land masses, and places the polar regions where we can see them. That way, our cruise ship and oil tanker captains can navigate more easily through the clear, blue Arctic Ocean, unimpeded by any tiresome ice-pack.
I particularly love that oil companies want to use the new Arctic Ocean sea lanes to transport their oil to market faster. Is it irony, or some form of rare, extinction-level stupidity that only comes around one every few thousand years? Hard to tell. But I digress.
The Standard Mercator Projection
As most of you will know, the standard world map in use on most walls today is an example of the spherical Mercator projection we are all familiar with. The problem with this projection is that near the poles it distorts any landmass, so that past 80º north or south of the equator it doesn’t really make much sense. You know, where Greenland looks half the size of Africa, when it’s actually only about a third as big as Australia. That’s because it’s an attempt to project a sphere onto a plane (OK, a vertical cylinder), using a straight line from the centre of the earth through every point on the earth’s surface, with the amount of distortion tending to infinity the closer you get to the north or south poles.
There are various assumptions in using this transform, among which is the fact that the earth is a perfect sphere, which it isn’t – it’s slightly oblate (squashed down) at the poles or, if you prefer, bulging at the equator, because of its spin. Putting these limitations aside, how do we relocate these points of maximum distortion to the open sea, away from any landmass, so that our polar zones can be seen more clearly? What if we moved our map poles to the equator, 90º E and W of Greenwich? This would place them in the Pacific Ocean, near the Galapagos Islands, and between Sri Lanka and Sumatra, in the Indian Ocean.
The Problem
This is nothing new, and it’s called the Transverse Mercator projection. But given a standard Mercator projection map, how do you create a Transverse Mercator projection? And how do you do it in Python?
Wikipedia gives us the equations. They are not trivial, and involve some pretty mind-bending trigonometry, but they are no match for Python’s Numpy package. There are various versions of the forward and reverse Transverse Mercator equations, but the ones I’ve used are shown below.
These are the forward equations, for calculating the Transverse Mercator map position of a given point of latitude and longitude:
And these are the inverse equations, for building the output map, line by line, from points on the normal Mercator input map:
The variables λ (lambda) and φ (phi) are the original latitude and longitude, while x and y are their Cartesian coordinates on your transformed map. The coefficients kºa are where you put all your normalising radian-to-pixel position calculations, depending on how large you want your output map to be. This includes conversion from normal Cartesian (x, y) coordinates to the (x, y) coordinates of an image, which go from top left to bottom right.
The forward equations are used to calculate how each point on a Mercator projection would be translated onto a new projection. This can be used to test the equations. But we also need the reverse equations to assemble the output image line by line, as if we’re working back from a resulting image that has yet not been created. If we only used the forward equations, there would be gaps in the output image. But by working back from an imaginary output image, it is as if we are asking “for each output coordinate (x,y) on the output image, which pixel from the input image would I use?” In other words, what (lat, long) coordinate does this transformed (x,y) coordinate represent?
Step One: Testing the Equations With Matplotlib
First, before writing any code to transform an actual map, the forward equations had to be tested. The goal was to see what the equations did, to find out what the range of output values was, and to check that the output map was not flipped or inverted in any way. This was achieved using Python’s Numpy package, combined with Python’s Matplotlib graphics plotting package. To see what the equations did, the program was first run on a set of latitude and longitudes for world cities, plus some islands and oceans, with the lines of latitude and longitude superimposed.
This was the result:
Any physicist or electrical engineer would notice that the transformed lines of latitude and longitude resemble the lines of force in an electromagnetic field. But for our purposes, the important thing is that even though the world’s cities have been translated from their normal positions, they are still in the same relative positions, albeit with the Northern Hemisphere on the left. From this we can discern how the transformed map would look, and what the equations are doing. We are now ready to test the equations on an actual map image.
The code to obtain the above plot is included below:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
This little program tests the forward transverse Mercator equations for
converting traditional latitude and longitude coords to their tranverse Mercator
equivalent. The result is a plot of the transformed positions of several cities,
islands and oceans.
Equations are from: en.wikipedia.org/wiki/
Transverse_Mercator_projection#Direct_transformation_formulae
Created on 5 May 2018
Author: matta_idlecoder at protonmail dot com
"""
import numpy as np
import matplotlib.pyplot as plt
from collections import namedtuple
import timeit
timing = True
def new_y_coord(lat, longtde):
"""Converts a lat, longitude location to a new transformed longitude
"""
lat_in_rads = lat / 180 * np.pi
long_in_rads = longtde / 180 * np.pi
new_y_rad = 0.5 * np.log((1 + (np.sin(long_in_rads) * np.cos(lat_in_rads)))
(1 - (np.sin(long_in_rads) * np.cos(lat_in_rads))))
new_y_deg = new_y_rad * 180 / np.pi
return new_y_deg
def new_x_coord(lat, longtde):
"""Converts a lat, longitude location to a new transformed latitude
"""
lat_in_rads = lat / 180 * np.pi
long_in_rads = longtde / 180 * np.pi
new_x_rad = -(np.arctan((np.tan(lat_in_rads)) / (np.cos(long_in_rads))))
new_x_deg = new_x_rad * 180 / np.pi
"""this code fixes the arctan problem of mapping only back to
-90 < arctan > 90. It does this by checking the position on the world map.
"""
if longtde < -90 or longtde > 90:
if lat > 0:
new_x_deg -= 180
elif lat < 0:
new_x_deg += 180
return new_x_deg
def plot_meridians():
"""Plots the transformed lines of latitude and longitude
"""
plt.figure(figsize=(15, 12), dpi=80)
plt.xlim(-180.0, 180.0)
plt.xticks(range(-90, 91, 90))
plt.ylim(-175.0, 175.0)
plt.yticks(range(-180, 181, 180))
plt.xlabel('Equator')
plt.ylabel('0 = International Dateline & \nGreenwich Meridian')
plt.title("Transverse Mercator Projection with Pole at 0'N, 90'W")
long_x, long_y = [], []
# plot the transformed lines of longitude:
for longitude in range(-180, 181, 10): # for ea vertical line of longitude
for latitude in range(-90, 91, 1): # print every point
if not latitude and abs(longitude) == 90: # don't div by 0 at poles
continue # don't calculate the point
x = round(new_x_coord(latitude, longitude), 2)
y = round(new_y_coord(latitude, longitude), 2)
long_x.append(x)
long_y.append(y)
plt.plot(long_x, long_y, 'g.', linewidth=0.5, label='Lines of Longitude')
Lat_x, Lat_y = [], []
# plot the transformed lines of latitude:
for latitude in range(-90, 91, 10): # for each horizontal line of latitude
for longitude in range(-180, 181, 1): # print every point
if not latitude and abs(longitude) == 90: # don't div by 0 at poles
continue # don't calculate the point
x = round(new_x_coord(latitude, longitude), 2)
y = round(new_y_coord(latitude, longitude), 2)
Lat_x.append(x)
Lat_y.append(y)
plt.plot(Lat_x, Lat_y, 'b--', linewidth=0.5, label='Lines of Latitude')
plt.legend(loc='upper right')
return
def plot_places():
"""Plots major world cities and places after a Shifted Mercator Transform
"""
City = namedtuple('City', 'name lat_long text_coords colour')
Ocean = namedtuple('Ocean', 'name lat_long text_coords')
Places = [
City(name='North Pole', lat_long=(90, 0), text_coords=(-130, 10),
colour='yellow'),
City(name='South Pole', lat_long=(-90, 0), text_coords=(110, -5),
colour='yellow'),
City(name='Greenwich meridian', lat_long=(5,0), text_coords=(-18, -15),
colour='green'),
City(name='Edinburgh', lat_long=(56, -4), text_coords=(-45, -5),
colour='blue'),
City(name='Sydney', lat_long=(-33, 152), text_coords=(130, 40),
colour='red'),
City(name='Rio', lat_long=(-22, -43), text_coords=(0, -40),
colour='black'),
City(name='Singapore', lat_long=(2, 104), text_coords=(-170, 90),
colour='black'),
City(name='S India', lat_long=(7, 77), text_coords=(-57, 97),
colour='orange'),
City(name='Quito', lat_long=(-2, -78), text_coords=(-20, -130),
colour='blue'),
City(name='Buenos Aires', lat_long=(-35, -56), text_coords=(40, -70),
colour='brown'),
City(name='Tierra del Fuego', lat_long=(-55, -65), text_coords=(100, -40),
colour='yellow'),
City(name='Caracas', lat_long=(10, -67), text_coords=(-50, -70),
colour='orange'),
City(name='Mexico City', lat_long=(19, -98), text_coords=(-145, -100),
colour='red'),
City(name='Los Angeles', lat_long=(33, -118), text_coords=(-145, -75),
colour='black'),
City(name='New York', lat_long=(42, -73), text_coords=(-55, -45),
colour='violet'),
City(name='Tokyo', lat_long=(36, 140), text_coords=(-160, 45),
colour='green'),
City(name='Auckland', lat_long=(-37, 175), text_coords=(150, -5),
colour='red'),
City(name='Moscow', lat_long=(56, 37), text_coords=(-52, 25),
colour='red'),
City(name='Cape Town', lat_long=(-34, 19), text_coords=(0, 10),
colour='brown'),
City(name='Cairo', lat_long=(30, 32), text_coords=(-25, 30),
colour='red'),
City(name='Baghdad', lat_long=(33, 45), text_coords=(-30, 50),
colour='blue'),
City(name='Bangkok', lat_long=(14, 101), text_coords=(-115, 120),
colour='yellow'),
City(name='Beijing', lat_long=(41, 115), text_coords=(-100, 60),
colour='orange'),
City(name='Kerguelen Is (Fr)', lat_long=(-49, 70), text_coords=(63, 63),
colour='black'),
City(name='SGeorgia', lat_long=(-55, -35), text_coords=(32, -11),
colour='yellow')]
Oceans = [
Ocean(name='Indian Ocean', lat_long=(-20, 80), text_coords=(40, 135)),
Ocean(name='N Pacific Ocean', lat_long=(30, -170),
text_coords=(-175, -38)),
Ocean(name='S Pacific Ocean', lat_long=(-30, 140),
text_coords=(-129, -80))]
for place in Places:
city_x, city_y = new_x_coord(*place.lat_long), \
new_y_coord(*place.lat_long)
plt.annotate(place.name, xy=(city_x, city_y), xytext=place.text_coords,
arrowprops=dict(facecolor=place.colour, shrink=0.05))
for watery_place in Oceans:
ocean_x, ocean_y = new_x_coord(*watery_place.lat_long), \
new_y_coord(*watery_place.lat_long)
plt.annotate(watery_place.name, xy=(ocean_x, ocean_y),
xytext=watery_place.text_coords)
return
if __name__ == '__main__':
if timing:
zerotime = timeit.default_timer()
plot_meridians()
plot_places()
plt.show()
if timing:
now = timeit.default_timer()
print('That took', round(now - zerotime, 2), 'seconds.\n')
Step Two: Converting Real Maps
Now that we know how the equations work, we can run the reverse equations on every pixel position in an imaginary output image, noting the pixel colour in the original input map for each output position. The plan was to create an intermediate 2D pixel position mapping, from which the output image would be assembled, pixel by pixel, line by line, from the original map image.
Using the original map at the top of the page, here is the result, now with our old North and South Poles lining up along our new equator:
As expected, the Northern Hemisphere is now on the left, and the Southern Hemisphere on the right. There is a lot of distortion of the Central and South American regions, and India looks a little bigger, but the polar regions are the areas of minimum distortion, which is what we wanted. We can also see that Canada, Scandinavia and Russian are much closer to their true sizes, now that they are nearer the ‘equator’.
Some Sample Outputs
The great thing about this program is that you can run it on any of your favourite Mercator maps, magically transforming them into transverse Mercator projections. Here is a sample of some other transformations I ran, on maps I found online:
How to Use This Program
The code to create the above transverse Mercator projections is given below. This can be used to convert any normal Mercator projection of the world to a transverse Mercator projection:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
{0}: redraws a given standard Mercator projection world map with the poles
moved to 90E & 90W at the equator. Will run on any image format or size of input
map. Assumes maps stop at 85º N/S.
Equations are explained at:
en.wikipedia.org/wiki/Transverse_Mercator_projection#Inverse_transformation_formulae
The following flags are optional:
-i input_map
The name of the input map filename. No extension required: program will
auto-detect the image type. With no -i flag, it will look for an
input file called WorldMap.png
-o output_map
The name you want to give your transformed output map. With no -o flag,
it will output to a file called ShiftedPoleMap.jpg
Requires an image type suffix, which can be any of the standard image
file types (.gif .bmp .jpg .jpeg .ppm .png .tif or .tiff).
Defaults to a .jpg output if none given.
-d
don't display maps when running the program
-h
Prints this __doc__ string
-t
Time the program execution
Example command line:
python3 {0} -i myMap.png -o myNewMap.png -t
Or, if you cd to the same folder and make this file executable:
./{0} -i myMap.png -o myNewMap.png -t
Created on 10 May 2018
Author: matta_idlecoder at protonmail dot com
"""
from PIL import Image
import numpy as np
import sys
import getopt
import timeit
import os
timing_it = True
displaying_images = True
def inv_merc_lat(translat_in_rads, translong_in_rads):
"""Converts a transverse mercator latitude to its original latitude.
"""
orig_lat_rad = -np.arcsin(
(np.sin(translong_in_rads)) / (np.cosh(translat_in_rads)))
return orig_lat_rad
def inv_merc_long(trans_lat_rads, trans_long_rads):
"""Converts a transverse mercator longitude to its original longitude.
"""
orig_long_rads = np.arctan(
np.sinh(trans_lat_rads) / np.cos(trans_long_rads))
if (trans_long_rads > np.pi / 2) or (trans_long_rads < -np.pi / 2): if trans_lat_rads >= 0:
orig_long_rads += np.pi
elif trans_lat_rads < 0:
orig_long_rads -= np.pi
return orig_long_rads
def inv_transform_pixel(x_pix_num, x_max, y_max, transformed_lat):
"""Calculates the new pixel position in the transformed image
"""
trans_longitude = round((float(x_pix_num) / x_max * 2 * np.pi) - np.pi, 2)
orig_latitude = round(inv_merc_lat(transformed_lat, trans_longitude), 3)
orig_longitude = round(inv_merc_long(transformed_lat, trans_longitude), 3)
original_x = int((orig_longitude + np.pi) / (2 * np.pi) * x_max)
original_y = int((orig_latitude - np.pi / 2) / (-np.pi) * y_max)
# fixes potential rounding and index errors:
original_x -= 1 if original_x == x_max else 0
original_y -= 1 if original_y == y_max else 0
return original_y, original_x
def inv_transform_line(map_line_num, x_max, y_max):
"""Returns a correlation between each pixel and its transformed equivalent
"""
# Not an image line, a mapping between pixels and their transformed positions:
line_pixel_mapping = np.zeros((x_max, 2),
dtype=np.uint16) # pix # will exceed 256
y_norm = float(map_line_num) / y_max
trans_lat = round(y_norm * (-35 * np.pi / 18) + (17 * np.pi / 18), 2)
for x_coord in range(x_max): # for each pixel on the line
line_pixel_mapping[x_coord] = inv_transform_pixel(x_coord, x_max, y_max,
trans_lat)
return line_pixel_mapping
def inv_transform_map(map_input_name, map_output_name):
"""Builds a new transformed map image line by line, pixel by pixel
"""
try:
input_image = Image.open(map_input_name)
X_PIC_SIZE, Y_PIC_SIZE = input_image.size
if displaying_images:
input_image.show()
except FileNotFoundError:
print("\nInput map {} isn't there. Check the filename is correct.\n".
format(map_input_name))
sys.exit(2)
print("\nRunning the map transform on a {}x{} map...".
format(X_PIC_SIZE, Y_PIC_SIZE))
line_pixel_mapping = np.zeros((Y_PIC_SIZE, X_PIC_SIZE, 2), dtype=np.uint16)
for map_image_line in range(Y_PIC_SIZE):
line_pixel_mapping[map_image_line] = inv_transform_line(map_image_line,
X_PIC_SIZE,
Y_PIC_SIZE)
# Create transformed image data:
output_image_data = np.zeros((Y_PIC_SIZE, X_PIC_SIZE, 3), dtype=np.uint8)
for line in range(Y_PIC_SIZE):
for pixel in range(X_PIC_SIZE):
new_pix_line, new_pix_x = line_pixel_mapping[line][pixel]
# The [:3] is to read only the RGB input image data if it's RGBA:
output_image_data[line][pixel] = input_image.getpixel(
(int(new_pix_x), int(new_pix_line)))[:3]
# Create transformed image:
output_image = Image.fromarray(output_image_data, 'RGB')
if displaying_images:
output_image.show()
if not map_output_name:
map_output_name = "ShiftedPoleMap.png"
output_image.save(map_output_name)
return
def usage(prog_name):
"""Prints a short command line guide, using whatever filename this prog has.
"""
print("""\nUsage: from the command line type:
{} [-i input_filename] [-o output_filename.ext] [-dht]
If you don't want to use flags, give your input file the name WorldMap.png
and your output filename will be given the name ShiftedPoleMap.jpg
Use the -t flag to time your map transformation.
Use the -h flag for a detailed flag listing.
""".format(prog_name))
return
def main(argv):
global displaying_images, timing_it
input_map_name = 'WorldMap.png'
output_map_name = 'ShiftedPoleMap.jpg'
try:
opts, args = getopt.getopt(argv, "i:o:dht")
except getopt.GetoptError as err:
print("\n", str(err))
usage(sys.argv[0])
sys.exit(2)
for opt, arg in opts:
if opt == '-d':
displaying_images = False
elif opt == '-i':
input_map_name = arg
elif opt == '-o':
output_map_name = arg
ext = os.path.splitext(output_map_name)[-1].lower()
if ext not in ".gif .bmp .jpg .jpeg .ppm .png .tif .tiff".split():
output_map_name += '.jpg' # smallest file format
print('Unknown output file extension. Creating a JPEG file.')
elif opt == '-h':
print("\n", __doc__.format(sys.argv[0][2:]))
sys.exit(2)
elif opt == '-t':
timing_it = True
else:
assert False, "unhandled option"
if timing_it:
start = timeit.default_timer()
inv_transform_map(input_map_name, output_map_name)
if timing_it:
stop = timeit.default_timer()
print('\nFinished the map transform in', round((stop - start), 2),
'seconds.\n')
return
if __name__ == '__main__':
main(sys.argv[1:])
I’ve written the code to be launched from the command line to convert any standard Mercator projection map. It can be launched in a number of ways. The simplest way, if you don’t like using command line flags, is to rename your input map WorldMap.jpg, to copy it to the same folder as the code, and then run the code as follows:
$ python3 MercatorInvTrans.py
Alternatively, you can make your file executable and run it as follows:
$ chmod +x MercatorInvTrans.py $ ./MercatorInvTrans.py
If you have some different maps you want to run it on, you can feed the program their names, using the standard Unix command line format, which will work in OS/X, Linux or Windows. The -i flag indicates that the following filename is the input file, and -o the output filename you want to use:
$ ./MercatorInvTrans.py -i myMap.png -o myNewMap.jpg
The program will save to any image format, which you select on the command line by giving your output file your preferred extension. if you omit the input image extension, it will auto-detect the image type. Make sure you install Python’s pillow package to get the latest version of the Python Image Library (PIL), which will handle all the image formatting. I’ve tested it with .bmp, .gif, .jpeg, .jpg, .png, .ppm. .tif and and .tiff extensions, and they all work. If you leave out the output file extension it will default to a .jpg file, which is usually the smallest image data format.
When you run the program, the input and output files are displayed automatically. If you don’t want them shown, use the -d for ‘don’t display the files’.
The other flag I’ve provided is -t, for timing the transformation. The calculation time for a small image of 720×360 pixels was a few seconds. For bigger images of 10k x 5k pixels or larger, it might take a few minutes. You can use the -h flag at any time to get these instructions.
Those unfamiliar with command lines should know that you can place the flags in any order on the command line, and that the Python getopt() utility in main() will work it out, e.g:
$ ./MercatorInvTrans.py -t -o myNewMap.png —i myMap
One final tip is to make sure your input Mercator map is symmetrical north and south about the equator, to the same maximum latitudes. The program assumes the input maps cut off at +/-85º. Note that many maps show up to 85º north, but only 75º south. If this is the case, you should use your favourite image editor to add some more Antarctica to the bottom of your input image, before you run the transform.
Further Reading
The Transverse Mercator transform gives a great excuse to explore the features of Matplotlib and the Python Image Library (PIL). But for those of you interested in exploring the map projection side of things further, you should know that there are some great map transform resources available online.
One great site that compares how different projections work is the SnipView site. Another is the site put together by Drew Roos called Mercator:Extreme. For map nerds, it’s an awesome resource. Written in Javascript, it allows you to drag your map’s north pole around anywhere on the planet, to see how the Mercator projection changes. As you play with this you realise you will always have to distort a landmass, wherever you place the poles – it’s just a matter of choosing which region to magnify.
I’ll leave you with a vision of the near future:
some truly interesting info , well written and broadly speaking user genial.
Cheers, Lon.
An excellent post, congratulations !!
Heya i’m for the first time here. I came across this board and I find It really useful & it helped
me out much. I hope to give something back and aid others like
you helped me.