python_scripts/add_trees_to_open_street_map/.gitignore

-cache
-*.html
+cache
\ No newline at end of file
-*.csv
-*.png
\ No newline at end of file

python_scripts/add_trees_to_open_street_map/add_trees.py

@@ -17,6 +17,8 @@ import overpy
 from pyproj import Transformer
 from shapely import LineString, geometry, wkt
 from shapely.ops import transform
+import pandas as pd
+import geopandas as gpd
 
 from tree import Forest
 from import_existing_trees import get_existing_forest
@@ -31,6 +33,7 @@ from import_existing_trees import get_existing_forest
 WKT = "POLYGON((9.947021 49.803063, 9.947011 49.800917, 9.955025 49.800810, 9.955110 49.803019, 9.947021 49.803063))"
 # Replace with None if no existing tree should be imported
 EXISTING_TREES = 'existing_trees/Trees_ideal_2_20240227.shp'
+# EXISTING_TREES = 'existing_trees/baumkataster/Baum.shp'
 # Fellbach
 # WKT = "POLYGON((9.271353 48.811327, 9.271911 48.809010, 9.272147 48.807187, 9.275838 48.807173, 9.275602 48.806749, 9.276138 48.806325, 9.277683 48.806424, 9.277319 48.812514, 9.275581 48.811991, 9.271353 48.811327))"
 EPSG_ID = 25832
@@ -46,11 +49,13 @@ IGNORE_ROADS = set(['primary', 'unclassified', 'secondary',
 
 
 SCRIPT_DIR = Path(__file__).resolve().parent
+OUTPUT_DIR = SCRIPT_DIR / 'output'
+Bounds = namedtuple("Bounds", "W S E N")
 
 
 def load_region(wkt_polygon):
     region = wkt.loads(wkt_polygon)
-    bounds = namedtuple("Bounds", "W S E N")(*region.bounds)
+    bounds = Bounds(*region.bounds)
     return region, bounds
 
 
@@ -139,13 +144,14 @@ def place_trees(forest, ways, region, to_local, tree_distance, min_distance_2) -
                                             color='#DFFF00',
                                             type='Fake Tree',
                                             description='Tilia tomentosa',
-                                            radius=3
+                                            diameter=6,
+                                            source='add_trees.py'
                                             )
 
     return forest
 
 
-def plot_trees(bounds, forest, tree_distance) -> None:
+def plot_trees(bounds: Bounds, forest: Forest, tree_distance: float) -> None:
     print("Exporting diagram...")
     tree_xs, tree_ys, colors = forest.xs_ys_cs
     plt.scatter(tree_xs, tree_ys, s=2, c=colors)
@@ -154,7 +160,8 @@ def plot_trees(bounds, forest, tree_distance) -> None:
     plt.title(f"{bounds}\nTree distance : {tree_distance} m")
     plt.gcf().set_size_inches(15, 10)
     plt.savefig(
-        SCRIPT_DIR / f"{get_basename(bounds)}.png", bbox_inches='tight', dpi=300)
+        OUTPUT_DIR / f"{get_basename(bounds)}.png", bbox_inches='tight', dpi=300)
+    print("  DONE!")
 
 
 def export_map(bounds, forest, epsg_id) -> None:
@@ -167,7 +174,7 @@ def export_map(bounds, forest, epsg_id) -> None:
         lon, lat = to_wgs84.transform(tree.x, tree.y)
         folium.Circle(
             location=[lat, lon],
-            radius=tree.radius or 2,  # [m], when defined
+            radius=tree.radius,  # [m]
             color="black",
             weight=1,
             fill_opacity=0.9,
@@ -186,21 +193,48 @@ def export_map(bounds, forest, epsg_id) -> None:
         control=True
     ).add_to(interactive_map)
 
-    interactive_map.save(f"{get_basename(bounds)}_trees.html")
+    interactive_map.save(OUTPUT_DIR / f"{get_basename(bounds)}_trees.html")
     print("  DONE!")
 
 
 def export_csv(bounds, forest, wkt_polygon, tree_distance, min_distance, epsg_id) -> None:
     print("Exporting CSV...")
-    with open(SCRIPT_DIR / f"{get_basename(bounds)}_trees.csv", "w") as csv:
+    with open(OUTPUT_DIR / f"{get_basename(bounds)}_trees.csv", "w") as csv:
         csv.write(f"# Fake trees for; {wkt_polygon}\n")
         csv.write(f"# Tree distance along roads; {tree_distance}; [m]\n")
         csv.write(f"# Minimum allowed distance between trees; {min_distance}; [m]\n")
         csv.write(f"# EPSG; {epsg_id}\n")
-        csv.write("# X; Y; Type; Description; Radius\n")
-        csv.write("# [m]; [m]; [?]; [?]; [m]\n")
+        csv.write("# X; Y; Type; Description; Radius; Source\n")
+        csv.write("# [m]; [m]; [-]; [-]; [m]; [-]\n")
         for tree in forest:
-            csv.write(f"{tree.x};{tree.y};{tree.type};{tree.description};{tree.radius}\n")
+            csv.write(f"{tree.x};{tree.y};{tree.type};{tree.description};{tree.radius};{tree.source}\n")
+
+    print("  DONE!")
+
+
+def export_shapefile(bounds: Bounds, forest: Forest, tree_distance: float, epsg_id: str) -> None:
+    print("Exporting shapefile")
+
+    data = [{
+        'x': t.x, 'y': t.y,
+        'Bezeichnun': t.description, 'Baumart': t.type,
+        'Baumhöhe': t.height, 'Kronenbrei': t.diameter,
+        'Quelle': t.source,
+        }
+        for t in forest]
+
+    df = pd.DataFrame.from_dict(data)
+    gdf = gpd.GeoDataFrame(
+        df.drop(columns=['x', 'y']),
+        geometry=gpd.points_from_xy(df.x, df.y), crs=epsg_id
+    )
+
+    print(gdf)
+
+    basename = get_basename(bounds)
+    shp_dir = OUTPUT_DIR / basename
+    shp_dir.mkdir(exist_ok=True)
+    gdf.to_file(shp_dir / f"trees.shp")
 
     print("  DONE!")
 
@@ -224,8 +258,8 @@ def main(wkt_polygon, epsg_id, tree_distance, min_distance, import_tree_shp) ->
 
     plot_trees(bounds, forest, tree_distance)
     export_map(bounds, forest, epsg_id)
-    export_csv(bounds, forest, wkt_polygon,
-               tree_distance, min_distance, epsg_id)
+    export_csv(bounds, forest, wkt_polygon, tree_distance, min_distance, epsg_id)
+    export_shapefile(bounds, forest, tree_distance, epsg_id)
 
 
 if __name__ == "__main__":

python_scripts/add_trees_to_open_street_map/existing_trees/baumkataster/Baum.cpg

+UTF-8
\ No newline at end of file

python_scripts/add_trees_to_open_street_map/existing_trees/baumkataster/Baum.prj

+PROJCS["ETRS89.UTM-32N",GEOGCS["LL-ETRF89",DATUM["ETRF89",SPHEROID["GRS1980",6378137.000,298.25722210],TOWGS84[0.0000,0.0000,0.0000,0.000000,0.000000,0.000000,0.00000000]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["false_easting",500000.000],PARAMETER["false_northing",0.000],PARAMETER["central_meridian",9.00000000000000],PARAMETER["scale_factor",0.9996],PARAMETER["latitude_of_origin",0.000],UNIT["Meter",1.00000000000000]]
\ No newline at end of file

python_scripts/add_trees_to_open_street_map/import_existing_trees.py

+from pathlib import Path
 import geopandas as gpd
 from tree import Tree, Forest
 
@@ -9,8 +10,10 @@ def get_existing_forest(shp_input):
         point = tree_row.geometry
         trees.append(Tree(point.x, point.y,
                           description=tree_row.Bezeichnun,
-                          radius=tree_row.Kroneradi,
-                          type=tree_row.Baumart
+                          diameter=tree_row.Kronenbrei,
+                          type=tree_row.Baumart,
+                          trunk_diameter=tree_row.Stammumfan,
+                          source=Path(shp_input).name
                           ))
     return Forest(trees)
 

python_scripts/add_trees_to_open_street_map/output/.gitignore

+*.png
+*.html
+*.csv
+*.cpg
+*.dbf
+*.shp
+*.prj
+*.shx
+*.geojson
\ No newline at end of file

python_scripts/add_trees_to_open_street_map/tree.py

@@ -7,11 +7,13 @@ import kdtree
 class Tree:
     x: float
     y: float
+    diameter: float
     z: float = 0
     height: float = None
-    radius: float = None
+    trunk_diameter: float = None
     type: str = None
     description: str = '?'
+    source: str = '?'
     color: str = 'green'
 
     def __len__(self):
@@ -20,8 +22,12 @@ class Tree:
     def __getitem__(self, i):
         return [self.x, self.y][i]
 
+    @property
+    def radius(self):
+        return self.diameter / 2
+
     def __str__(self):
-        return f"{self.type} ({self.description}), {self.radius or '?'} m (X={self.x:.1f}, Y={self.y:.1f})"
+        return f"{self.type} ({self.description}), {self.radius} m (X={self.x:.1f}, Y={self.y:.1f})"
 
 
 class Forest(UserList):