Frolian's blog standing on the shoulders of PyGiants

• About
• Archive

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

Thu, 30 Oct 2014

Geography 3D Printing

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.

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).

from IPython.display import Image
Image(url="http://thingiverse-production.s3.amazonaws.com/renders/d7/de/ac/dc/f4/sakurajima_preview_featured.jpg")

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.

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):
    elvtn_args.update({
    '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']:
        elevationArray.append(resultset['elevation'])

    return elevationArray

# 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:

from pylab import *
%matplotlib inline

plot(elevationArray)

[<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.

elvtn_args = {}
elvtn_args.update({
    '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.

response.keys()

['status', 'results']

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

response['status']

'OK'

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

response['results'][:5]

[{'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.

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)
xlabel('latitude')
ylabel('longitude')
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.

top_left = (48.960092, 2.158060) #latitude, longitude
bottom_right = (48.802953, 2.456751)

from mpl_toolkits.basemap import Basemap

# setup Lambert Conformal basemap.
figure(figsize=(10, 10))
m = Basemap(width=1000000, height=1000000, projection='lcc',
            resolution='l',
            lat_0=47.,lon_0=2.)
# draw coastlines
m.drawcoastlines()
# 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.
m.drawmapboundary(fill_color='aqua')
# fill continents, set lake color same as ocean color.
m.fillcontinents(color='coral',lake_color='aqua')
m.drawcountries()
# 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)
m.drawmeridians(meridians,labels=[True,False,False,True])


# 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],
         arrowprops=dict(arrowstyle="->",
                connectionstyle="angle3,angleA=0,angleB=-90"),
         fontsize=15) 

<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.

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.

import time

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)
        time.sleep(0.5)
        response = simplejson.load(urllib.urlopen(url))
    data.append(response)

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.

def build_elevationArray():
    elevationArray = []
    for result_set in data:
        for result_point in result_set['results']:
            elevationArray.append(result_point['elevation'])
    return array(elevationArray).reshape((200, 200))

elevationArray = build_elevationArray()

Y, X = meshgrid(latitudes, longitudes) # needed for nicely displaying the contours

figure(figsize=(15, 15))
imshow(elevationArray, cmap='cubehelix', extent=(min_longitude, max_longitude, latitudes.min(), latitudes.max()))
colorbar(shrink=0.5)
title('terrain elevation (m)')
xlabel('longitude')
ylabel('latitude')

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

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

figure(figsize=(15, 7))
contour(X, Y, elevationArray.T, 20, cmap='cubehelix')
colorbar()
title('terrain elevation (m)')
xlabel('longitude')
ylabel('latitude')

<matplotlib.text.Text at 0x8509450>

Building a 3D model from the height data

In the following cells, we will build a 3D model from the height data. To achieve this, we will extract the geometrical data from the contour sets we just used in the previous plot and extrude them in the $z$ direction. The extrusion will be done by generating code for the OpenSCAD modeller. This will allow us to directly convert the model to the meshed STL format used by a 3D printer.

The reason I am generating code for OpenSCAD is that it's really easy to use. For instance, extruding a triangle of height and width equal to 100 mm in the $z$ direction over 10 mm can be expressed by the following line:

linear_extrude(height = 10) polygon(points=[[0,0],[100,0],[0,100]]);

Image(filename='files/OpenSCAD.PNG')

First, we generate the contour we want to work with.

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

Next, we extract a collection of paths for each isosurface (whose existence I have discovered by reading this link - before I thought that I should use "segments"...). Each path is made of vertices, which we can plot. The lines obtained by this will be extruded in our 3D model.

def get_vertices(CS, isosurface_index):
    lc = CS.collections[isosurface_index]
    vertices = []
    for path in lc.get_paths():
        vertices.append(path.vertices)
    return vertices

These vertices can be plotted:

def plot_isosurface(isosurface_index):
    figure(1, figsize=(10, 4))
    
    subplot(121)
    title("height {0} meters".format(CS.levels[isosurface_index]))
    vertices = get_vertices(CS, isosurface_index)
    for v in vertices:
        plot(v[:, 0], v[:, 1])
    xlim(longitudes.min(), longitudes.max())
    ylim(latitudes.min(), latitudes.max())
    
    subplot(122)
    contour(X, Y, elevationArray.T, CS.levels, cmap='cubehelix')

plot_isosurface(2)

This can even be made interactive.

from IPython.html.widgets import interact

interact(plot_isosurface,
         isosurface_index=(0, len(CS.collections) - 1))

Now that we have all the vertices we need in our mesh, we can build up a 3D extension of those. We need to convert the coordinates from each point of the segment from latitude longitude coordinates to some coordinates appropriate to the OpenSCAD modelling.

We compute the extent of our model below.

lon_mini = longitudes.min()
lat_mini = latitudes.min()
dx = longitudes.max() - longitudes.min()
dy = latitudes.max() - latitudes.min()

dx, dy

(0.29869100000000026, 0.15713900000000081)

Next, we define a scaling factor that will permit to compute the final dimensions of our object. Our scaling factor will be the target $dx$ (longitude extent) of our object, expressed in millimeters.

scaling_factor = 100. #mm

This allows us to write some scaling functions.

scale_lon = lambda lon: scaling_factor / dx * (lon - lon_mini) 
scale_lat = lambda lat: scaling_factor / dx * (lat - lat_mini)

From this, we can define a function that scales the coordinates of a point and returns some OpenSCAD compatible string:

array2str = lambda arr: "[{0:.5f}, {1:.5f}]".format(scale_lon(arr[0]), scale_lat(arr[1]))

Using the previous function, we can define a function that loops over all points and returns an extruded version using the polygon command of OpenSCAD. The function below takes into account a starting elevation due to the fact that we will stack extrusions on top of each other to build the final 3D model.

def vertices2extrusion(vertices, start_height, extrusion_height):
    polygon = "[" + ", ".join([array2str(point) for point in vertices]) + "]"
    return "translate([0, 0, {0}]) linear_extrude(height={1}) polygon({2});".format(start_height, extrusion_height, polygon)

vertices2extrusion(get_vertices(CS, 2)[0], 0, 10)

'translate([0, 0, 0]) linear_extrude(height=10) polygon([[0.00000, 49.52006], [0.35937, 49.70117], [0.50251, 49.83891], [0.79369, 49.96554], [1.00503, 50.05929], [1.19791, 50.22991], [1.44856, 50.49427], [1.50754, 50.54141], [1.73315, 50.75864], [1.89860, 51.02301], [2.01005, 51.17476], [2.10872, 51.28738], [2.38678, 51.55175], [2.51256, 51.75968], [2.56659, 51.81611], [2.78789, 52.08048], [3.01508, 52.23854], [3.14294, 52.34485], [3.01508, 52.44844], [2.71004, 52.60922]]);'

vertices2extrusion(get_vertices(CS, 2)[0], 20, 10)

'translate([0, 0, 20]) linear_extrude(height=10) polygon([[0.00000, 49.52006], [0.35937, 49.70117], [0.50251, 49.83891], [0.79369, 49.96554], [1.00503, 50.05929], [1.19791, 50.22991], [1.44856, 50.49427], [1.50754, 50.54141], [1.73315, 50.75864], [1.89860, 51.02301], [2.01005, 51.17476], [2.10872, 51.28738], [2.38678, 51.55175], [2.51256, 51.75968], [2.56659, 51.81611], [2.78789, 52.08048], [3.01508, 52.23854], [3.14294, 52.34485], [3.01508, 52.44844], [2.71004, 52.60922]]);'

Now that we have our basic functions, we can apply this to all segments in our contours and output this to a 3D file. We will use define a z scaling factor much in the same way that we scaled our latitude and longitude above.

scaling_factor_z = 50. / 200. # target size in z is 50 mm
scale_z = lambda z: scaling_factor_z * z 

def generate_model_code(CS, filename='files/model.scad'):
    """ generates OpenSCAD model code and writes it to file.
    """
    model = ""
    for index in range(len(CS.collections)): # loop over every isosurface
        current_height = CS.levels[index]
        if index != (len(CS.collections) - 1):
            current_extrusion = CS.levels[index + 1] - CS.levels[index]
        else:
            current_extrusion = CS.levels[index] - CS.levels[index - 1]
        vertices_list = get_vertices(CS, index)
        
        # loop for one layer
        layer = ""
        for v in vertices_list:
            layer += vertices2extrusion(v, 
                                        scale_z(current_height),
                                        scale_z(current_extrusion))
        model += layer
    
    # write to file
    with open(filename, 'w') as f:
        f.write(model)
                                    

generate_model_code(CS)

At this point, we have a working model. However, I noticed that we don't have a flat plate upon which the whole structure rests. So at this point I'm going back one step and manually setting the borders of our data to zero so that the isocontours we get are naturally bounded.

elevationArray[:, 0] = elevationArray[:, -1] = 0
elevationArray[0, :] = elevationArray[-1, :] = 0
CS = contour(X, Y, elevationArray.T, 20)

As one can see, there now are borders on the edges of the model. This naturally builds a support plate for our 3D model.

generate_model_code(CS, filename='files/model_bounded.scad')

Another experiment we can do is to spatially filter the data with a Gaussian filter. This will smooth all contours.

from scipy.ndimage import gaussian_filter

elevationArray = build_elevationArray() # rebuild data
elevationArray = gaussian_filter(elevationArray, 2) # filter it

elevationArray[:, 0] = elevationArray[:, -1] = 0 # make sure there's a nice border
elevationArray[0, :] = elevationArray[-1, :] = 0

CS = contour(X, Y, elevationArray.T, 50) # build contour set

generate_model_code(CS, filename='files/model_bounded_smoothed.scad') # write model

Finally, we can display a little gallery of the generated models we now have.

from IPython.display import display
for f in ['files/OpenSCAD_model.PNG', 'files/OpenSCAD_model_bounded.PNG', 'files/OpenSCAD_model_bounded_smoothed.PNG']:
    display(Image(filename=f))

Final and last step: print it using a 3D printer!

Conclusions

In this post, we have learned how to use the Google elevation API to retrieve height data for a grid of geographical coordinates centered about Paris. This data was used in conjunction with the powerful numeric Python tools to produce a smooth mesh that was rendered using the OpenSCAD language. This 3D model can be exported to the STL format, ready for 3D printing.

The tools developed within this post are generic and could thus make it possible to produce a similarly rendered 3D model of any region of the world.

This post was entirely written using the IPython notebook. You can see a static view or download this notebook with the help of nbviewer at 20141027_ParisTopology.ipynb.

---

This article has been tagged with the following terms:

Geography 3D Printing

---

• 
• 
• 
• 
• 

---

Comments