From 315494ea33fe39875afe8c5cd27dbc75db3fea1e Mon Sep 17 00:00:00 2001 From: Brent Schroeter Date: Sun, 18 May 2025 21:54:59 -0700 Subject: [PATCH] add option for water cutouts --- README.md | 1 + src/__init__.py | 110 ++++++++++++++++++++++++++++++++++++------------ 2 files changed, 84 insertions(+), 27 deletions(-) diff --git a/README.md b/README.md index afb8b57..abac21b 100644 --- a/README.md +++ b/README.md @@ -34,6 +34,7 @@ docker run --rm -it -v "$(pwd)/data:/app/data" geo3dp:latest \ Other available options include: - `--square`: crops the model area to a square centered within the bounding box - `--baseplate-h`: adjusts the height of the rectangular base +- `--water-cutouts`: creates through holes where there are water features ## Parallel Rendering diff --git a/src/__init__.py b/src/__init__.py index 2d18160..d350d6c 100644 --- a/src/__init__.py +++ b/src/__init__.py @@ -100,7 +100,7 @@ def get_building_height(building): """ Infer the height, in meters, of a building from OSM building data. """ - DEFAULT_METERS_PER_LEVEL = 4 + DEFAULT_METERS_PER_LEVEL = 3.5 DEFAULT_LEVELS_PER_BUILDING = 3 if building["height"] is not None: return float(building["height"]) @@ -109,6 +109,16 @@ def get_building_height(building): return DEFAULT_METERS_PER_LEVEL * DEFAULT_LEVELS_PER_BUILDING +def remove_holes(shape): + """ + Remove holes from a shapely Polygon or MultiPolygon. Always returns a + MultiPolygon. + """ + if shape.geom_type == "Polygon": + shape = MultiPolygon([shape]) + return MultiPolygon([Polygon(geom.exterior) for geom in shape.geoms]) + + class ElevationMap: def __init__(self, tiff_file_path): self._gdal = gdal.Open(tiff_file_path) @@ -149,6 +159,7 @@ class Renderer: osm_path: str elev_grid: np.array baseplate_h: float + cut_out_water: bool output: str def render_chunk(self): @@ -187,9 +198,46 @@ class Renderer: origin=(0, 0), ) - model = solid2.cube([chunk_bounds[2] - chunk_bounds[0], chunk_bounds[3] - chunk_bounds[1], baseplate_h]) + model = solid2.cube([ + chunk_bounds[2] - chunk_bounds[0], + chunk_bounds[3] - chunk_bounds[1], + baseplate_h, + ]) ROAD_HEIGHT = 0.25 + + print(f"[{self.output}] 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(chunk_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)) + # Render roads on slopes by over-extruding them and then calculating their # intersection with a slightly offset elevation grid drive_mask = None @@ -219,6 +267,28 @@ class Renderer: else: drive_mask += drive_mask_cube + if self.cut_out_water: + print(f"[{self.output}] Constructing water features...") + natural = osm.get_natural() + water_bodies = natural[(natural["natural"] == "water") | natural["water"].notna()] + water_polygons = [] + for _, body in water_bodies.iterrows(): + shape = convert_shape_from_world_coords(body["geometry"]) + shape = shape.simplify(0.1) + shape = remove_holes(shape) + shape = shape.intersection(chunk_rect) + if shape.area > 2: + water_polygons.append(remove_holes(shape)) + water_polygons = union_all(water_polygons) + # Buffering slightly helps to remove artifacts at the tile edge + water_polygons = water_polygons.buffer(0.05, cap_style="flat") + if water_polygons.geom_type == "Polygon": + water_polygons = MultiPolygon([water_polygons]) + for shape in water_polygons.geoms: + model -= solid2.linear_extrude(100)( + to_scad_polygon(shape), + ) + print(f"[{self.output}] Constructing networks...") drive_polygons = [] for shape in drive_net["geometry"]: @@ -229,31 +299,9 @@ class Renderer: if drive_polygons.geom_type == "Polygon": drive_polygons = MultiPolygon([drive_polygons]) for shape in drive_polygons.geoms: - model += solid2.linear_extrude(100)(to_scad_polygon(shape)).translate(0, 0, baseplate_h) * drive_mask - - print(f"[{self.output}] 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(chunk_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)) + model += solid2.linear_extrude(100)( + to_scad_polygon(shape), + ) * drive_mask print(f"[{self.output}] Rendering...") render_stl(model, self.output) @@ -304,6 +352,7 @@ if __name__ == "__main__": parser.add_argument( "--square", action="store_true", + default=False, help="Constrains the model to a square with dimension --edge-len", ) parser.add_argument( @@ -312,6 +361,12 @@ if __name__ == "__main__": default=5.0, help="Depth of the base plate in millimeters", ) + parser.add_argument( + "--water-cutouts", + action="store_true", + default=False, + help="Creates holes for water features", + ) parser.add_argument( "--output", help="File path for STL output", @@ -410,6 +465,7 @@ if __name__ == "__main__": int(elev_grid.shape[1]/args.slices*j):int(elev_grid.shape[1]/args.slices*(j+1)), ], baseplate_h=args.baseplate_h, + cut_out_water=args.water_cutouts, output=f"/tmp/chunk.{i}.{j}.stl", )) with Pool(cpu_count()) as pool: