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