import argparse from math import asin, cos, pi, sin, sqrt import numpy as np import solid2 from osgeo import gdal from pyrosm import OSM from shapely import affinity, Polygon, union_all def haversine_dist(lat1, lng1, lat2, lng2): """ Calculate distance between two globe coordinates in meters. """ EARTH_RADIUS = 6378137 # in meters lat1_rads = lat1 * pi / 180 lat2_rads = lat2 * pi / 180 lng1_rads = lng1 * pi / 180 lng2_rads = lng2 * pi / 180 d_lat = lat2_rads - lat1_rads d_lng = lng2_rads - lng1_rads hav_theta = sin(d_lat / 2)**2 + cos(lat1_rads) * cos(lat2_rads) * sin(d_lng / 2)**2 return 2 * EARTH_RADIUS * asin(sqrt(hav_theta)) def to_scad_polygon(shapely_polygon): """ Generate a solid2 SCAD polygon object from a shapely Polygon. """ points = list(shapely_polygon.exterior.coords) paths = [[i for i, _ in enumerate(shapely_polygon.exterior.coords)]] for interior in shapely_polygon.interiors: paths.append([len(points) + i for i in range(len(interior.coords))]) points += list(interior.coords) return solid2.polygon(points, paths) def get_building_height(building): """ Infer the height, in meters, of a building from OSM building data. """ DEFAULT_METERS_PER_LEVEL = 5 DEFAULT_LEVELS_PER_BUILDING = 3 if building["height"] is not None: return float(building["height"]) if building["building:levels"] is not None: return float(building["building:levels"]) * DEFAULT_METERS_PER_LEVEL return DEFAULT_METERS_PER_LEVEL * DEFAULT_LEVELS_PER_BUILDING class ElevationMap: def __init__(self, tiff_file_path): self._gdal = gdal.Open(tiff_file_path) self._transform = self._gdal.GetGeoTransform() # Gets band at index 1 from the data self._band = self._gdal.GetRasterBand(1) def get_elevation(self, lat, lng): """ Get elevation, in meters above sea level, for a point within the loaded dataset. """ x = int((lng - self._transform[0]) / self._transform[1]) y = int((lat - self._transform[3]) / self._transform[5]) if 0 <= x < self._gdal.RasterXSize and 0 <= y < self._gdal.RasterYSize: return self._band.ReadAsArray(x, y, 1, 1)[0, 0] raise IndexError("coordinates are outside of elevation map area") if __name__ == "__main__": parser = argparse.ArgumentParser() parser.add_argument( "--osm", help=".osm.pbf file", required=True, ) parser.add_argument( "--topo", help=".tif file from https://portal.opentopography.org/raster?opentopoID=OTSRTM.082015.4326.1", required=True, ) parser.add_argument( "--bbox", help=",,,", required=True, ) parser.add_argument( "--edge-len", type=float, default=100.0, help="Long edge length of output STL bounding box", ) parser.add_argument( "--output", help="File path for STL output", required=True, ) parser.add_argument( "--baseplate-h", type=float, default=5.0, help="Depth of the base plate in millimeters", ) parser.add_argument( "--circular", action="store_true", ) args = parser.parse_args() elev_map = ElevationMap(args.topo) osm = OSM(args.osm) solid2.set_global_fn(48) drive_net = osm.get_network(network_type="driving") natural = osm.get_natural() buildings = osm.get_buildings() map_bounds = [float(coord) for coord in args.bbox.split(",")] map_width = map_bounds[2] - map_bounds[0] map_height = map_bounds[3] - map_bounds[1] # Figure out how much to squish the map's coordinate system to make it # square meters_per_unit_lat = haversine_dist( map_bounds[1], map_bounds[0], map_bounds[3], map_bounds[0], ) / map_height meters_per_unit_lng = haversine_dist( map_bounds[1], map_bounds[0], map_bounds[1], map_bounds[2], ) / map_width lat_lng_ratio = meters_per_unit_lat / meters_per_unit_lng map_width_meters = map_width * meters_per_unit_lng map_height_meters = map_height * meters_per_unit_lat aspect_ratio = map_width_meters / map_height_meters if aspect_ratio > 1: # Wide model_bounds = (0, 0, args.edge_len, args.edge_len / aspect_ratio) else: # Tall model_bounds = (0, 0, args.edge_len * aspect_ratio, args.edge_len) model_rect = Polygon([ (0, 0), (0, model_bounds[3]), (model_bounds[2], model_bounds[3]), (model_bounds[2], 0), ]) # Scale factor from real-world meters to the model coordinate system scale_factor = model_bounds[2] / map_width_meters print("Map bounds:", map_bounds) print("Meters per degree latitude:", meters_per_unit_lat) print("Meters per degree longitude:", meters_per_unit_lng) print("Aspect ratio:", aspect_ratio) print("Output bounds:", model_bounds) print( "Scale factor (millimeter in model per meter in world):", f"{scale_factor:.4g}", ) def convert_shape_from_world_coords(shape): # Move the coordinate system so that the origin is at the bottom left # of the map rectangle, then scale the shape in x and y directions return affinity.scale( affinity.translate(shape, -map_bounds[0], -map_bounds[1]), model_bounds[2] / map_width, model_bounds[3] / map_height, origin=(0, 0), ) model = solid2.cube([model_bounds[2] - model_bounds[0], model_bounds[3] - model_bounds[1], args.baseplate_h]) print("Calculating elevations...") # Slice the model area into a rectilinear grid to approximate topography ELEV_RES = 30 elev_grid = np.zeros((ELEV_RES, ELEV_RES)) for i in range(ELEV_RES): for j in range(ELEV_RES): elev_grid[i, j] = elev_map.get_elevation( map_bounds[1] + (map_bounds[3] - map_bounds[1]) / ELEV_RES * j, map_bounds[0] + (map_bounds[2] - map_bounds[0]) / ELEV_RES * i, ) min_world_elev = elev_grid.min() elev_grid = (elev_grid - min_world_elev) * scale_factor ROAD_HEIGHT = 0.25 # Render roads on slopes by over-extruding them and then calculating their # intersection with a slightly offset elevation grid drive_mask = None for i in range(ELEV_RES): for j in range(ELEV_RES): if elev_grid[i, j] > 0: model += solid2.cube([ model_bounds[2] / ELEV_RES, model_bounds[3] / ELEV_RES, elev_grid[i, j] + args.baseplate_h, ]).translate( model_bounds[2] / ELEV_RES * i, model_bounds[3] / ELEV_RES * j, 0, ) drive_mask_cube = solid2.cube([ model_bounds[2] / ELEV_RES, model_bounds[3] / ELEV_RES, elev_grid[i, j] + args.baseplate_h + ROAD_HEIGHT, ]).translate( model_bounds[2] / ELEV_RES * i, model_bounds[3] / ELEV_RES * j, 0, ) if drive_mask is None: drive_mask = drive_mask_cube else: drive_mask += drive_mask_cube print("Constructing networks...") drive_polygons = [] for shape in drive_net["geometry"]: shape = convert_shape_from_world_coords(shape).simplify(1).buffer(0.25) drive_polygons.append(shape.intersection(model_rect)) drive_polygons = [shape for shape in drive_polygons if not shape.is_empty] for shape in union_all(drive_polygons).geoms: model += solid2.linear_extrude(100)(to_scad_polygon(shape)).translate(0, 0, args.baseplate_h) * drive_mask print("Constructing buildings...") for i, building in buildings.iterrows(): shape = convert_shape_from_world_coords(building["geometry"]) shape = shape.buffer(0.125) # Simplifying building shapes after scaling can speed up rendering by # an order of magnitude shape = shape.simplify(0.05) shape = shape.intersection(model_rect) if not shape.is_empty: height = get_building_height(building) * scale_factor building_center = building["geometry"].centroid # Computing elevation from the elevation map can result in weird # artifacts, so calculate based on the less granular precomputed # grid instead elev = elev_grid[ min(max(int((building_center.x - map_bounds[0]) / map_width * ELEV_RES), 0), ELEV_RES - 1), min(max(int((building_center.y - map_bounds[1]) / map_height * ELEV_RES), 0), ELEV_RES - 1), ] if shape.geom_type == "Polygon": model += solid2.linear_extrude(args.baseplate_h + elev + ROAD_HEIGHT + height)(to_scad_polygon(shape)) else: for shape in shape.geoms: model += solid2.linear_extrude(args.baseplate_h + elev + ROAD_HEIGHT + height)(to_scad_polygon(shape)) if args.circular: radius = min(model_bounds[2], model_bounds[3]) / 2 model *= solid2.translate(model_bounds[2] / 2, model_bounds[3] / 2, 0)( solid2.cylinder(100, r=radius), ) print("Rendering... (grab a cup of coffee)") solid2.render_to_stl_file(model, args.output)