A 3D model of the Paris area using the Google elevation API and OpenSCAD

Update November 3, 2014: I printed the model using an Ultimaker. You can find it on Thingiverse. Here's what it looks like when printed. Paris Relief

In this post, we will use the elevation API provided by Google maps to obtain some height information about a section of Paris, and then build an 3D model to be printed from the data. The end result we are trying to achieve should look much the same as the image shown below (the source of the image is this Thingiverse thing).

In [1]:
from IPython.display import Image

Let's get started.

Getting the source height data

Sample data

Our first step is to get the height data we need for our 3D representation. We're going to use Google data that we can access using the elevation API. The elevation API usage is described here. Users of the free API can do the following:

  • 2,500 requests per 24 hour period.
  • 512 locations per request.
  • 25,000 total locations per 24 hour period.
  • 5 requests per second.

Below, we're executing the (sligthly modified) sample script found on the API documentation page to get some elevation data along a line defined by start and end coordinates.

In [2]:
import simplejson
import urllib

ELEVATION_BASE_URL = 'https://maps.googleapis.com/maps/api/elevation/json'
CHART_BASE_URL = 'http://chart.apis.google.com/chart'

def getElevation(path="36.578581,-118.291994|36.23998,-116.83171",samples="100", **elvtn_args):
    'path': path,
    'samples': samples})

    url = ELEVATION_BASE_URL + '?' + urllib.urlencode(elvtn_args)
    response = simplejson.load(urllib.urlopen(url))

    # Create a dictionary for each results[] object
    elevationArray = []

    for resultset in response['results']:

    return elevationArray
In [3]:
# Mt. Whitney
startStr = "36.578581,-118.291994"
# Death Valley
endStr = "36.23998,-116.83171"

pathStr = startStr + "|" + endStr

elevationArray = getElevation(pathStr)

To visualize the data we have just downloaded, we can plot it:

In [4]:
from pylab import *
%matplotlib inline
In [5]:
[<matplotlib.lines.Line2D at 0x583e0d0>]

This "line data" will be used in the next section to obtain all the data we need for a 3D model. Before we move on, I'd like to explore the detailed results obtained by the API request.

In [6]:
elvtn_args = {}
    'path': pathStr,
    'samples': "200"})
url = ELEVATION_BASE_URL + '?' + urllib.urlencode(elvtn_args)
response = simplejson.load(urllib.urlopen(url))

As one can see below, response is a dictionary containing the status and the results related to the request we just made.

In [7]:
['status', 'results']

The status reports whether the request worked or not. Here, it is:

In [8]:

The results are a list of dictionaries containing the coordinate points.

In [9]:
[{'elevation': 4350.96044921875,
  'location': {'lat': 36.578581, 'lng': -118.291994},
  'resolution': 305.4064636230469},
 {'elevation': 3765.516845703125,
  'location': {'lat': 36.57692408822881, 'lng': -118.2846242328207},
  'resolution': 305.4064636230469},
 {'elevation': 3697.862548828125,
  'location': {'lat': 36.57526672284261, 'lng': -118.27725478197},
  'resolution': 305.4064636230469},
 {'elevation': 3629.592529296875,
  'location': {'lat': 36.57360890388832, 'lng': -118.2698856475017},
  'resolution': 305.4064636230469},
 {'elevation': 3324.3681640625,
  'location': {'lat': 36.57195063141282, 'lng': -118.2625168294698},
  'resolution': 305.4064636230469}]

These coordinates can be extracted to produce a sample 3D plot.

In [10]:
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
xs = [item['location']['lat'] for item in response['results']]
ys = [item['location']['lng'] for item in response['results']]
zs = [item['elevation'] for item in response['results']]
ax.plot3D(xs, ys, zs)
ax.set_zlabel('elevation (m)')
<matplotlib.text.Text at 0x5a53bd0>

Defining a grid of elevations and retrieving heights

Now that we know how to retrieve elevation data from a given path defined by start and end points, we can extend this method to work with a grid. In the next cells, we will discretize a portion of the globe roughly corresponding to the location of Paris (although any portion of map could do) and gather its elevation data.

First, we define the portion of the map we want to use by its coordinates (in fact, its "top left" and "bottom right"). We also plot it using Basemap to make sure that our bounding coordinates are correctly selected.

In [11]:
top_left = (48.960092, 2.158060) #latitude, longitude
bottom_right = (48.802953, 2.456751)
In [12]:
from mpl_toolkits.basemap import Basemap
In [13]:
# setup Lambert Conformal basemap.
figure(figsize=(10, 10))
m = Basemap(width=1000000, height=1000000, projection='lcc',
# draw coastlines
# draw a boundary around the map, fill the background.
# this background will end up being the ocean color, since
# the continents will be drawn on top.
# fill continents, set lake color same as ocean color.
# meridians on bottom and left
parallels = np.arange(45, 55, 5.)
# labels = [left,right,top,bottom]
m.drawparallels(parallels, labels=[False,True,True,False])
meridians = np.arange(-2.5, 10, 2.5)

# computing remaining points of bounding rectangle
top_right = (top_left[0], bottom_right[1])
bottom_left = (bottom_right[0], top_left[1])

for start, end in zip([top_left, top_right, bottom_right, bottom_left],
                      [top_right, bottom_right, bottom_left, top_left]):
    x0, y0 = m(start[1], start[0])
    x1, y1 = m(end[1], end[0])
    m.plot([x0, x1], [y0, y1], '-bo')

region_center = (0.5 * (top_left[0] + bottom_right[0]),
                 0.5 * (top_left[1] + bottom_right[1]))
xy_region_center = m(region_center[1], region_center[0])

annotate("region to extract", [xy_region_center[0], xy_region_center[1]], 
         [xy_region_center[0] + 50000, xy_region_center[1] + 50000],
<matplotlib.text.Annotation at 0x65075b0>

This looks like what we expect: a rectangle of data located in the northern part of France.

The next step is to discretize a grid based on the coordinates that we have been using so far so that we can retrieve elevation data for each grid point to create a height map. We create a grid of latitudes and longitudes based on the bounding box coordinates.

In [14]:
min_longitude, max_longitude = top_left[1], bottom_right[1]
latitudes = linspace(top_left[0], bottom_right[0], 200)
longitudes = linspace(min_longitude, max_longitude, 200)

And now we loop over all trajectories extracting 200 points on each line using the Google elevation API.

In [15]:
import time
In [16]:
start_t = time.time()
data = []
for ind, lat in enumerate(latitudes):
    path = "{0},{1}|{2},{3}".format(lat, min_longitude, lat, max_longitude)
    elvtn_args = {'path': path,
                  'samples': '200'}
    url = ELEVATION_BASE_URL + '?' + urllib.urlencode(elvtn_args)
    response = simplejson.load(urllib.urlopen(url))
    while not response['status'] == 'OK':
        print "response was not ok for request(index: {0}, path: {1}), will retry shortly.".format(ind, path)
        response = simplejson.load(urllib.urlopen(url))

That's it. We have now downloaded the data we need for our 3D modelling task.

Plotting the height data

Below, we loop over the results in our dataset and append them to a list. We then reshape the data into a matrix for easier handling. Finally, we build a coordinate mesh and plot the elevation data and some contour surfaces.

In [17]:
def build_elevationArray():
    elevationArray = []
    for result_set in data:
        for result_point in result_set['results']:
    return array(elevationArray).reshape((200, 200))
In [18]:
elevationArray = build_elevationArray()
In [19]:
Y, X = meshgrid(latitudes, longitudes) # needed for nicely displaying the contours
In [20]:
figure(figsize=(15, 15))
imshow(elevationArray, cmap='cubehelix', extent=(min_longitude, max_longitude, latitudes.min(), latitudes.max()))
title('terrain elevation (m)')

CS = contour(X, Y, elevationArray.T, 20)

This looks quite beautiful. We can also just plot the contours of this map.

In [21]:
figure(figsize=(15, 7))
contour(X, Y, elevationArray.T, 20, cmap='cubehelix')
title('terrain elevation (m)')
<matplotlib.text.Text at 0x8509450>