My 7 Hills: GeoJSON to STL

tl;dr: We reconstruct an STL mesh from GeoJSON data by extracting the GeoJSON shape coordinates into a point cloud and applying Poisson reconstruction using Open3d.

Deeply unsatisfied with my first attempt at printing GeoJSON hills, I came up with a three-step process that produces high-fidelity reproductions using only open source software. This process has three steps:

1. Slice the large geojson map using Geopandas.
2. Create an XYZ point cloud of the geojson geometry (leveraging Shapely and Geopandas).
3. Build a TriangleMesh for the XYZ point cloud using Open3D.
4. Extrude a solid from this mesh using Blender.

For this iteration, we’ll try to model the cliffside including Telegraph hill. If you want to try this yourself, here’s my script and the SF GeoJSON data.

Telegraph Hill is an iconic cliff-sided elevation topped by Coit Tower.

Step 1: Slicing the earth

GeoPandas extends Pandas dataframes with geodata and spacial primitives. It allows us to clip shapes to our Longitude, Latitude designated window of interest, convert coordinates to UTM X and Y meters, and more.

First, let’s slice the SF geojson map to only include features (and aspects of features) within the coordinates of Telegraph Hill:

#!/usr/bin/env python3
import geopandas
from shapely.geometry import Polygon

outname = "telegraph_hill"
y1, x1 = (37.8053307084282, -122.40853235179131)
y2, x2 = (37.79953896610929, -122.40128101638189)

gdf = geopandas.read_file('geojson/sf.geojson')
polygon = Polygon([(x1, y1), (x1, y2), (x2, y2), (x2, y1), (x1, y1)])
gdf_clipped = gdf.clip(polygon)
gdf_clipped.to_file(f"geojson/{outname}.geojson", driver='GeoJSON')
gdf_clipped = gdf_clipped.explode(index_parts=True)  #Line Segments only!

Step 2: Generating a Point Cloud from GeoJSON

Step 1 gave us a “GeoDataFrame” of shapes within our desired coordinates. We now need to convert these relevant shapes to a point cloud, and then into a triangle mesh for 3d printing.

We will continue to use geopandas to translate the coordinates from Long, Lat to X, Y meters from origin, and join these with the Z data from our source data’s ‘elevation’ property. These datapoints then initialize an Open3d point cloud.

import open3d as open3d

CONV_FT_TO_M = 0.3048  # SFData provides elevation in feet :(

def compute_gdf_pointcloud(gdf: geopandas.GeoDataFrame, z_scale: float) -> o3d.geometry.PointCloud:
    """Compute the PointCloud of the GeoDataFrame."""
    min_x, min_y = (99999999,99999999)
    min_elevation = min([e for e in gdf['elevation']])
    
    for f in gdf['geometry']:
        for c in f.coords:
            if c[0] < min_x:
                min_x = c[0]
            if c[1] < min_y:
                min_y = c[1]

    logging.debug(f"min_x={min_x} min_y={min_y} min_elevation={min_elevation}")
    gdf['flat_coordinates'] = gdf['geometry'].combine(
        gdf['elevation'],
        (lambda g, e: [(float(c[0] - min_x),
                        float(c[1] - min_y),
                        float(e) * CONV_FT_TO_M * z_scale) for c in g.coords]))
    # Add these coordinates to an Open3d point cloud.
    pcd = o3d.geometry.PointCloud()
    pcd.points = o3d.utility.Vector3dVector()
    for coordinate_list in gdf['flat_coordinates']:
        pcd.points.extend(o3d.utility.Vector3dVector(coordinate_list))
    # Estimate normals.
    pcd.estimate_normals(search_param=o3d.geometry.KDTreeSearchParamHybrid(radius=0.01, max_nn=30))
    pcd.orient_normals_to_align_with_direction()
    return pcd

Step 3: Build a TriangleMesh for the XYZ point cloud

We now attempt Poisson mesh reconstruction from our point cloud with normals. The Poisson mesh reconstruction also returns density data for each inferred triangle, which allows us to trim extraneously inferred areas if necessary. (It also the only interesting coloration we have)

import open3d as o3d
import numpy as np
import matplotlib.pyplot as plt

def compute_poisson_mesh(pcd: o3d.geometry.PointCloud, depth: int) 
   -> o3d.geometry.TriangleMesh:
    """Compute the mesh of the point cloud.
    depth:    The depth of the octree used for the surface reconstruction determines resolution of the resulting triangle mesh.
    """
    poisson_mesh, densities = o3d.geometry.TriangleMesh.create_from_point_cloud_poisson(
        pcd, depth=depth, width=0, scale=1.1, linear_fit=True)
    # Color the supporting point density.
    densities = np.asarray(densities)
    graph_rate = (densities - densities.min()) / (densities.max() - densities.min())
    density_colors = plt.get_cmap('plasma')(graph_rate)
    density_colors = density_colors[:, :3]
    poisson_mesh.vertex_colors = o3d.utility.Vector3dVector(density_colors)
    # Trim any excess surface from the Poisson reconstruction.
    bbox = pcd.get_axis_aligned_bounding_box()
    p_mesh_crop = poisson_mesh.crop(bbox)
    p_mesh_crop.compute_triangle_normals()
    return p_mesh_crop

Altogether not bad! Inferred from the point density, the color of the mesh indicates how many points (from our pointcloud) indicate confidence in the given polygon. Since our process benefits from a rectangular mesh (see the next step), we don’t trim any excess area.

Step 4: Fixing and printing

We need to fix a few holes in our geometry and optionally create an efficient solid (mountains should be hollow). To fix holes, import your STL into Blender and select the vertices surrounding the hole in edit mode and hit “alt-F” to fill the hole with triangles.

A generalized solution for STL hole patching is left as an exercise for the reader.

We should now be able to extrude our mesh and use the Bisect tool to preserve the upper part. Check for non-manifold faces along the way, and recompute normals before exporting.

This is much nicer than what I was initially expecting! I hope this helps you in your journey, and happy hacking.