commit 316cbe065adf103244044b1fe31bfb4187c764af Author: Brent Schroeter Date: Sat May 17 00:14:02 2025 -0700 initial commit diff --git a/Dockerfile b/Dockerfile new file mode 100644 index 0000000..962b515 --- /dev/null +++ b/Dockerfile @@ -0,0 +1,16 @@ +FROM python:3.13.3-bookworm + +WORKDIR /app + +RUN apt-get update && apt-get install -y openscad libgdal-dev +RUN python3 -m pip install setuptools wheel markupsafe numpy + +COPY requirements.txt . + +RUN python3 -m pip install -r requirements.txt + +COPY src src + +VOLUME /app/data + +ENTRYPOINT ["python3", "src/__init__.py"] diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..6d706bf --- /dev/null +++ b/requirements.txt @@ -0,0 +1,5 @@ +geopandas==1.0.1 +gdal[numpy]==3.6.2 +pyrosm==0.6.2 +shapely==2.1.0 +solidpython2==2.1.1 diff --git a/src/__init__.py b/src/__init__.py new file mode 100644 index 0000000..73c7e02 --- /dev/null +++ b/src/__init__.py @@ -0,0 +1,264 @@ +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)