Using Python To Make A Map With Accurate Longitude/Latitude Axes From The Database At OpenStreetMap.Org

Introduction

In order to work with the map database at openstreetmap.org the website has extensive information and one can find very many software packages available for developers of different kinds of uses for maps. This author had a particular application for maps where the longitude and latitude grid is to be registered accurately with the map image. Furthermore, in developing the application it was desired to use a general purpose programming tool to accomplish this and to avoid downloading and learning to use any of the software packages listed at openstreetmap.org. In this case that programming tool was python3 and its various powerful libraries.

In order to obtain and work with the database one needs to understand Slippy maps. This web page describes everything to know about Slippy map tiles and how to get them.

Slippy Maps

A brief tutorial is needed but without the detail found in the referenced web page.

The openstreetmap.org database is constructed with tiles, each tile being a 256 pixel x 256 pixel .png image file. There are 20 directories where each directory contains all of the tiles that cover the world, where the world is a Mercator projection that one commonly sees when looking at any map.

The directories are called zoom levels, and are enumerated 0 thru 19. Directory zoom level 0 contains a single tile. Zoom level 1 contains 4 tiles, creating a 2 x 2 grid that covers the world. Greater resolution is found as the zoom level increases. Zoom level 12 consists of a 4096 x 4096 grid of tiles covering the world. The following image taken from the referenced web page depicts the grid for zoom level 2.

Tiles are enumerated left to right (x) and top to bottom (y). Note this is done in the image. The zoom level and x, y enumeration allows one to specify each particular png file.

Web clients can download tiles from the openstreetmap.org web servers by sending http requests. The format is: https://tile.openstreetmap.org/zoom/x/y.png

There are three subdomain servers from which one can also download tiles. In the url “tile” is replaced with “a.tile” or “b.tile” or “c.tile”. At this point it is important to note that openstreetmap.org has a policy to limit overloading its servers which occurs by specifying high zoom levels together with building a large area map. A map is constructed by assembling the tiles in a grid and clearly the zoom level defines the resolution of the map. However, the user should trade off resolution and make do with fewer tiles to avoid needing to download hundreds of tiles.

Lastly, it is most important to know how to translate a longitude/latitude coordinate into the particular tile designation in which that coordinate is contained, that is, to find zoom/x/y.png. For python code this calculation is as follows:

import math
def deg2num(lat_deg, lon_deg, zoom):
  lat_rad = math.radians(lat_deg)
  n = 2.0 ** zoom
  xtile = int((lon_deg + 180.0) / 360.0 * n)
  ytile = int((1.0 - math.asinh(math.tan(lat_rad)) / math.pi) / 2.0 * n)
  return (xtile, ytile)

Now to emphasize the point that is the reason for this post. All one knows after downloading a tile is that the longitude/latitude point specified in the above function lies somewhere in the obtained tile, but not where. Therefore, one can build a map with a grid of connected downloaded tiles around some point on earth, but to calibrate that map more work will be needed.

Building A Map With Python

A very good piece of python code derived from “BerndGit” at stackoverflow.com is used to show how to build a map in python.

#import matplotlib.pyplot as plt
import math
import requests
from io import BytesIO
from PIL import Image
import matplotlib
import sys
    

def deg2num(lat_deg, lon_deg, zoom):
  lat_rad = math.radians(lat_deg)
  n = 2.0 ** zoom
  xtile = int((lon_deg + 180.0) / 360.0 * n)
  ytile = int((1.0 - math.log(math.tan(lat_rad) + (1 / math.cos(lat_rad))) / math.pi) / 2.0 * n)
  return (xtile, ytile)
    
def num2deg(xtile, ytile, zoom):
  n = 2.0 ** zoom
  lon_deg = xtile / n * 360.0 - 180.0
  lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * ytile / n)))
  lat_deg = math.degrees(lat_rad)
  return (lat_deg, lon_deg)
   
def getImageCluster(lat_deg, lon_deg, delta_lat,  delta_long, zoom):
    headers = {"User-Agent":"Mozilla/5.0 (Macintosh; Intel Mac OS X 10_15_4) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/83.0.4103.97 Safari/537.36"}
    smurl = r"http://a.tile.openstreetmap.org/{0}/{1}/{2}.png"
    xmin, ymax =deg2num(lat_deg, lon_deg, zoom)
    xmax, ymin =deg2num(lat_deg + delta_lat, lon_deg + delta_long, zoom)
   
    Cluster = Image.new('RGB',((xmax-xmin+1)*256-1,(ymax-ymin+1)*256-1) ) 
    for xtile in range(xmin, xmax+1):
        for ytile in range(ymin,  ymax+1):
            try:
                imgurl = smurl.format(zoom, xtile, ytile)
                print("Opening: " + imgurl)
                imgstr = requests.get(imgurl, headers=headers)
                tile = Image.open(BytesIO(imgstr.content))
                Cluster.paste(tile, box = ((xtile-xmin)*256 ,  (ytile-ymin)*255))
            except: 
                print("Couldn't download image")
                tile = None
    Cluster.save("my_image2.png")  
    return Cluster
        a = getImageCluster(lat_deg, lon_deg, delta_lat,  delta_lon, zoom)
    
    
if __name__ == '__main__':

    top, bot = 37.843634, 37.620166 #finished image top latitude and bottom latitude      
    lat_deg=bot   
    delta_lat=top-bot
    lef, rgt = -121.963417, -121.553999 #finished image left longitude and right longitude  
    lon_deg=lef  
    delta_lon=rgt-lef    
    zoom = 12
    a = getImageCluster(lat_deg, lon_deg, delta_lat,  delta_lon, zoom)
   

A few items need to be explained. The function “num2deg” is not used here but will be shown what it does later. The library “matplotlib” is not used here but will be later. The inputs to this function are the latitudes and longitudes that define the boundaries of the desired map, that is, “top, bot” and “lef, rgt”, and also the desired zoom level. Imports from python libraries “io” and “PIL” do the heavy lifting. After the downloaded tiles are pasted together the final map is saved to a png file.

The last item needing explanation is the json “headers” key that is passed to the “requests.get” function. The need for this key is found at this web site. The following image shows the map image built with this code.

Now given how this works and explained above the bounding box, “top, bot” and “lef, rgt”, lies within this map image with its four corners near the respective four corners of the map image within 255 pixels of the x and y directions. Stated another way, this map image of pasted-together tiles is somewhat larger than the requested map, which is contained within.

So, how does one now apply a set of longitude and latitude axes to this created map? There are two ways and the following starts with the first way.

Finding the Longitude/Latitude Corners

At the same web post at stackoverflow.com another user “etna” builds upon the above code and uses the function “num2deg” to find the longitude and latitude coordinates of the four corners of map image that was created. That is, the map image that is somewhat larger than the image bounded by “top, bot” and “lef, rgt”. The following code shows how this is done.

import matplotlib.pyplot as plt
import math
import requests
from io import BytesIO
from PIL import Image
import matplotlib
import sys

def deg2num(lat_deg, lon_deg, zoom):
  lat_rad = math.radians(lat_deg)
  n = 2.0 ** zoom
  xtile = int((lon_deg + 180.0) / 360.0 * n)
  ytile = int((1.0 - math.log(math.tan(lat_rad) + (1 / math.cos(lat_rad))) / math.pi) / 2.0 * n)
  return (xtile, ytile)
    
def num2deg(xtile, ytile, zoom):
  """
  http://wiki.openstreetmap.org/wiki/Slippy_map_tilenames
  This returns the NW-corner of the square. 
  Use the function with xtile+1 and/or ytile+1 to get the other corners. 
  With x
"""
  n = 2.0 ** zoom
  lon_deg = xtile / n * 360.0 - 180.0
  lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * ytile / n)))
  lat_deg = math.degrees(lat_rad)
  return (lat_deg, lon_deg)
   
def getImageCluster(lat_deg, lon_deg, delta_lat,  delta_long, zoom):
    headers = {"User-Agent":"Mozilla/5.0 (Macintosh; Intel Mac OS X 10_15_4) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/83.0.4103.97 Safari/537.36"}
    smurl = r"http://a.tile.openstreetmap.org/{0}/{1}/{2}.png"
    xmin, ymax =deg2num(lat_deg, lon_deg, zoom)
    xmax, ymin =deg2num(lat_deg + delta_lat, lon_deg + delta_long, zoom)

    #the four corners of the map are calculated
    bbox_ul = num2deg(xmin, ymin, zoom) #bbox_ul[0]=top, bbox[1]=lef
    bbox_ll = num2deg(xmin, ymax + 1, zoom) #bbox_ll[0]=bot, bbox_ll[1]=lef
    print (bbox_ul, bbox_ll) #the left corners

    bbox_ur = num2deg(xmax + 1, ymin, zoom) #bbox_ur[0]=top, bbox_ur[1]=rgt
    bbox_lr = num2deg(xmax + 1, ymax +1, zoom) #bbox_ll[0]=bot, bbox_lr[1]=rgt
    print (bbox_ur, bbox_lr) #the right corners

   
    Cluster = Image.new('RGB',((xmax-xmin+1)*256-1,(ymax-ymin+1)*256-1) ) 
    for xtile in range(xmin, xmax+1):
        for ytile in range(ymin,  ymax+1):
            try:
                imgurl = smurl.format(zoom, xtile, ytile)
                print("Opening: " + imgurl)
                imgstr = requests.get(imgurl, headers=headers)
                tile = Image.open(BytesIO(imgstr.content))
                Cluster.paste(tile, box = ((xtile-xmin)*256 ,  (ytile-ymin)*255))
            except: 
                print("Couldn't download image")
                tile = None
    #Cluster.save("my_image2.png")  
    return Cluster, [bbox_ll[1], bbox_ll[0], bbox_ur[1], bbox_ur[0]] #the new lef, bot, rgt, top

if __name__ == '__main__':
    
    top, bot = 37.843634, 37.620166 #desired image top latitude and bottom latitude      
    lat_deg=bot   
    delta_lat=top-bot
    lef, rgt = -121.963417, -121.553999 #desired image left longitude and right longitude  
    lon_deg=lef  
    delta_lon=rgt-lef    
    zoom = 12
    a, bbox = getImageCluster(lat_deg, lon_deg, delta_lat,  delta_lon, zoom)     
    print ("input coordinates: " + str(top) + ", " + str(bot) + ", " + str(lef) + ", " + str(rgt))
    #save the original requested map
    lat = [ top, bot, top, bot ]
    lon = [ lef, lef, rgt, rgt ]

    #new bounding coordinates of downloaded, created map
    top, bot = bbox[3], bbox[1]  #latitudes
    lef, rgt = bbox[0], bbox[2]  #longitudes
    print ("resulting output coordinates: " + str(top) + ", " + str(bot) + ", " + str(lef) + ", " + str(rgt))

    img = a #this is now the base image to plot
    fig, ax = plt.subplots(figsize=(12,12))
    ax.scatter(lon, lat, alpha=0.5, c='b', s=10)
    ax.set_ylim(bot, top)
    ax.set_xlim(lef, rgt)
    #draw bounding lines showing the original map request
    plt.axvline(x = lon[0], color = 'r',  lw=0.9, ls='--' )
    plt.axvline(x = lon[2], color = 'r',  lw=0.9, ls='--' )
    plt.axhline(y = lat[0], color = 'r',  lw=0.9, ls='--' )
    plt.axhline(y = lat[1], color = 'r',  lw=0.9, ls='--' )
    ax.imshow(img, extent=(lef, rgt, bot, top), aspect="equal")
    filesave = "visual.png"
    plt.savefig(filesave, dpi=300, bbox_inches='tight', pad_inches=0)

In the function “num2deg” a comment has been added to describe what this function returns. In the function “getImageCluster” the calculation for the longitude/latitude coordinates for the four corners of the map is done. This produces new values for “top, bot” and “lef, rgt” which are then returned to the calling program along with the image “a”.

The code then proceeds to use matplotlib to develop a .png image of the map with accurate longitude/latitude axes. Additionally, the code shows the bounding box of the requested input values. The resulting image clearly shows that the requested box lies within the map created by pasting together the downloaded tiles even though the right edge is very near the map’s’s right edge in this particular case of values chosen for “top, bot” and “lef, rgt”.

Cropping the Map to Produce the Requested Map

There may be reason for the user to require a map precisely as requested by “top, bot” and “lef, rgt”, the longitude and latitude values. The way to do this is to crop the map image that results from the above code.

In the above code the python PIL library was used to create the image of the map “Cluster”. One now wants to crop this image using “Cluster.crop(….)”. To do this one needs to convert longitude values and latitude values to number of pixels since the crop function requires the pixel coordinates of the top left corner of the target image and its bottom right corner referenced to top left corner of “Cluster” as an origin.

To begin one calculates the pixel coordinates of the top left corner of “Cluster”. This corner is the top left tile that was downloaded, and in the above code its tile number coordinates are “xmin” and “ymin”.

#xmin and ymin calculated here in getImageCluster(...) 

xmin, ymax =deg2num(lat_deg, lon_deg, zoom)
xmax, ymin =deg2num(lat_deg + delta_lat, lon_deg + delta_long, zoom)

#in the example xmin=660 and ymin=1582
#now x,y can be calculated and these pixel coordinates become the origin for cropping

x, y = xmin * TILE_SIZE, ymin * TILE_SIZE #num pixels to left edge and num pixels to top edge

Now as stated earlier the “crop” function needs the pixel coordinates of the top left corner and the bottom right corner of the desired image. The longitude/latitude coordinates of these two points are the input request, “top, lef” and “bot, rgt”, respectively. To find the pixel coordinates a function similar to “deg2num” is needed and one can call it “deg2pix”.

#note that this function is almost identical to deg2num

TILE_SIZE = 256
def deg2pix(lat_deg, lon_deg, zoom):
  r = math.pow(2, zoom) * TILE_SIZE #total number of pixels at this zoom level
  lat_rad = math.radians(lat_deg)
  xt = int((lon_deg + 180.0) / 360.0 * r)
  y = int((1.0 - math.log(math.tan(lat_rad) + (1 / math.cos(lat_rad))) / math.pi) / 2.0 * r)
  return (x, y) #number of pixels to the point lon_deg and lat_deg, respectively

Now one can add to the function “getImageCluster(…)” and obtain the cropped image.

#x0, y0 corresponds to the top, lef point
#x1, y1 corresponds to the bot, rgt point

x0, y0 = deg2pix(lat_deg + delta_lat, lon_deg   , zoom)
x1, y1 = deg2pix( lat_deg, lon_deg + delta_long  , zoom)

Cluster = Cluster.crop((
        x0 - x,  # left
        y0 - y,  # top
        x1 - x,  # right
        y1 - y)) # bottom

In the numerical example this cropping operation is depicting in the next image.

Final Code-Cropped Image

Puttting the above together into final code and showing the resulting map:

import matplotlib.pyplot as plt
import numpy as np
import math
import requests
from io import BytesIO
from PIL import Image
import matplotlib
import sys

TILE_SIZE = 256
def deg2pix(lat_deg, lon_deg, zoom):
  r = math.pow(2, zoom) * TILE_SIZE
  lat_rad = math.radians(lat_deg)
  x = int((lon_deg + 180.0) / 360.0 * r)
  y = int((1.0 - math.log(math.tan(lat_rad) + (1 / math.cos(lat_rad))) / math.pi) / 2.0 * r)
  return (x, y)


def deg2num(lat_deg, lon_deg, zoom):
  lat_rad = math.radians(lat_deg)
  n = 2.0 ** zoom
  xtile = int((lon_deg + 180.0) / 360.0 * n)
  ytile = int((1.0 - math.log(math.tan(lat_rad) + (1 / math.cos(lat_rad))) / math.pi) / 2.0 * n)
  return (xtile, ytile)
    
def num2deg(xtile, ytile, zoom):
  n = 2.0 ** zoom
  lon_deg = xtile / n * 360.0 - 180.0
  lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * ytile / n)))
  lat_deg = math.degrees(lat_rad)
  return (lat_deg, lon_deg)
   
def getImageCluster(lat_deg, lon_deg, delta_lat,  delta_long, zoom):
    headers = {"User-Agent":"Mozilla/5.0 (Macintosh; Intel Mac OS X 10_15_4) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/83.0.4103.97 Safari/537.36"}
    smurl = r"http://a.tile.openstreetmap.org/{0}/{1}/{2}.png"
    xmin, ymax =deg2num(lat_deg, lon_deg, zoom) #tile coordinates of bottom left tile
    xmax, ymin =deg2num(lat_deg + delta_lat, lon_deg + delta_long, zoom)  #tile coordinates of top right tile
    #these are tile coordinates, not number of pixels
#    print ("xmin, ymax: " + str(xmin) + "," + str(ymax) + " tile coordinates of bottom left tile")
#    print ("xmax, ymin: " + str(xmax) + ","  + str(ymin) + " tile coordinates of top right tile")

    Cluster = Image.new('RGB',((xmax-xmin+1)*256-1,(ymax-ymin+1)*256-1) ) #the arguments are number of pixels 
    for xtile in range(xmin, xmax+1):
        for ytile in range(ymin,  ymax+1):
            try:
                imgurl = smurl.format(zoom, xtile, ytile)
                print("Opening: " + imgurl)
                imgstr = requests.get(imgurl, headers=headers)
                tile = Image.open(BytesIO(imgstr.content))
                Cluster.paste(tile, box = ((xtile-xmin)*256 ,  (ytile-ymin)*255))
            except: 
                print("Couldn't download image")
                tile = None

    x, y = xmin * TILE_SIZE, ymin * TILE_SIZE #num pixels to left edge and num pixels to top edge
    #cannot obtain x0, y0  and x1, y1 without defining deg2pix function whereas the deg2num function
    #has TILESIZE multiplied times the lat_rad  
    x0, y0 = deg2pix(lat_deg + delta_lat, lon_deg   , zoom)
    x1, y1 = deg2pix( lat_deg, lon_deg + delta_long  , zoom)
#    print ("x0, y0: " + str(x0) + "," +str(y0) + " number of pixels x and y to the top left corner of finished image") 
#    print ("x1, y1: " + str(x1) + "," +str(y1) + " number of pixels x and y to the bottom right corner of finished image") 
#    print ("x, y: " + str(x) + "," + str(y) + " actual x, y pixels of the left edge and top edge of the top left tile") 

    Cluster = Cluster.crop((
        x0 - x,  # left
        y0 - y,  # top
        x1 - x,  # right
        y1 - y)) # bottom
    return Cluster
    
    
    
if __name__ == '__main__':
    top, bot = 37.843634, 37.620166 #finished image top latitude and bottom latitude      
    lat_deg=bot   
    delta_lat=top-bot
    lef, rgt = -121.963417, -121.553999 #finished image left longitude and right longitude  
    lon_deg=lef  
    delta_lon=rgt-lef    
    zoom = 12
    a = getImageCluster(lat_deg, lon_deg, delta_lat,  delta_lon, zoom)
    img = a #this is now the base image to plot
    fig, ax = plt.subplots(figsize=(12,12))
    ax.set_ylim(bot, top)
    ax.set_xlim(lef, rgt)
    ax.imshow(img, extent=(lef, rgt, bot, top), aspect="equal")
    filesave = "visual_cropped.png"
    plt.savefig(filesave, dpi=300, bbox_inches='tight', pad_inches=0)