add option for water cutouts

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
 

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
+            model += solid2.linear_extrude(100)(
-
+                to_scad_polygon(shape),
-        print(f"[{self.output}] Constructing buildings...")
+            ) * drive_mask
-        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))
 
         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: