|
62 | 62 | plot_trials(circle_area_suite)
|
63 | 63 |
|
64 | 64 |
|
65 |
| -# ## Vancouver watershed benchmarks |
66 |
| -#= |
67 |
| -
|
68 |
| -Vancouver Island has ~1,300 watersheds. LibGEOS uses this exact data |
69 |
| -in their tests, so we'll use it in ours as well! |
70 |
| -
|
71 |
| -We'll start by loading the data, and then we'll use it to benchmark various operations. |
72 |
| -
|
73 |
| -=# |
74 |
| - |
75 |
| -# The CRS for this file is EPSG:3005, or as a PROJ string, |
76 |
| -# `"+proj=aea +lat_1=50 +lat_2=58.5 +lat_0=45 +lon_0=-126 +x_0=1000000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"` |
77 |
| -# TODO: this doesn't work with WellKnownGeometry. Why? |
78 |
| -wkt_geoms = LG.readgeom.(eachline("/Users/anshul/Downloads/watersheds.wkt"), (LG.WKTReader(LG.get_global_context()),)) |
79 |
| -vancouver_polygons = GI.getgeom.(wkt_geoms, 1); #.|> GO.tuples; |
80 |
| - |
81 |
| -import SortTileRecursiveTree as STR |
82 |
| -tree = STR.STRtree(vancouver_polygons) |
83 |
| -query_result = STR.query(tree, GI.extent(vancouver_polygons[1])) |
84 |
| - |
85 |
| -GO.intersects.((vancouver_polygons[1],), vancouver_polygons[query_result]) |
86 |
| - |
87 |
| -go_vp = GO.tuples.(vancouver_polygons[1:2]) |
88 |
| -@be GO.union($(go_vp[1]), $(go_vp[2]); target = $GI.PolygonTrait()) |
89 |
| -@be LG.union($(vancouver_polygons[1]), $(vancouver_polygons[2])) |
90 |
| - |
91 |
| -all_intersected = falses(length(vancouver_polygons)) |
92 |
| -accumulator = deepcopy(vancouver_polygons[1]) |
93 |
| -all_intersected[1] = true |
94 |
| -i = 1 |
95 |
| -# query_result = STR.query(tree, GI.extent(accumulator)) |
96 |
| -# for idx in query_result |
97 |
| -# println("Assessing $idx") |
98 |
| -# if !all_intersected[idx] && LG.intersects(vancouver_polygons[idx], accumulator) |
99 |
| -# println("Assimilating $idx") |
100 |
| -# result = LG.union(vancouver_polygons[idx], accumulator#=; target = GI.PolygonTrait()=#) |
101 |
| -# # @show length(result) |
102 |
| -# accumulator = result#[1] |
103 |
| -# all_intersected[idx] = true |
104 |
| -# end |
105 |
| -# end |
106 |
| -display(poly(vancouver_polygons[all_intersected]; color = rand(RGBf, sum(all_intersected)))) |
107 |
| -display(poly(accumulator)) |
108 |
| -@time while !(all(all_intersected)) && i < length(vancouver_polygons) |
109 |
| - println("Began iteration $i") |
110 |
| - query_result = STR.query(tree, GI.extent(accumulator)) |
111 |
| - for idx in query_result |
112 |
| - if !(all_intersected[idx] || !(LG.intersects(vancouver_polygons[idx], accumulator))) |
113 |
| - println("Assimilating $idx") |
114 |
| - result = LG.union(vancouver_polygons[idx], accumulator#=; target = GI.PolygonTrait()=#) |
115 |
| - # @show length(result) |
116 |
| - accumulator = result#[1] |
117 |
| - all_intersected[idx] = true |
118 |
| - end |
119 |
| - end |
120 |
| - display(poly(vancouver_polygons[all_intersected]; color = rand(RGBf, sum(all_intersected)))) |
121 |
| - println("Finished iteration $i") |
122 |
| - # wait_for_key("Press any key to continue to the next iteration.") |
123 |
| - i += 1 |
124 |
| -end |
125 |
| - |
126 |
| -fig = Figure() |
127 |
| -ax = Axis(fig[1, 1]; title = "STRTree query for polygons", aspect = DataAspect()) |
128 |
| -for (i, result_index) in enumerate(query_result) |
129 |
| - poly!(ax, vancouver_polygons[result_index]; color = Makie.wong_colors()[i], label = "$result_index") |
130 |
| -end |
131 |
| -Legend(fig[1, 2], ax) |
132 |
| -fig |
133 |
| - |
134 |
| - |
135 |
| - |
136 |
| -# TODO: |
137 |
| -# - Vancouver watersheds: |
138 |
| -# - Intersection on pre-buffered geometry |
139 |
| -# - Polygon union by reduction (perhaps pre-sort by border order, so we don't end up with useless polygons) |
140 |
| -# - Queries using STRTree.jl |
141 |
| -# - Potentially using a prepared geometry based approach to multithreaded reductive joining |
142 |
| -# - Implement multipolygon joining. How? Query intersection or touching for each individual geometry, |
143 |
| -# and implement a |
144 |
| - |
145 |
| - |
146 | 65 | ## Segmentization
|
147 | 66 |
|
148 | 67 |
|
@@ -180,9 +99,9 @@ n_points_values = round.(Int, exp10.(LinRange(0.7, 6, 15)))
|
180 | 99 | circle_union_suite["LibGEOS"][n_points] = @be LG.union($lg_circle_left, $lg_circle_right)
|
181 | 100 | end
|
182 | 101 |
|
183 |
| -plot_trials(circle_difference_suite) |
184 |
| -plot_trials(circle_intersection_suite) |
185 |
| -plot_trials(circle_union_suite) |
| 102 | +plot_trials(circle_difference_suite; legend_position = (2, 1), legend_kws = (; orientation = :horizontal, nbanks = 2)) |
| 103 | +plot_trials(circle_intersection_suite; legend_position = (2, 1), legend_kws = (; orientation = :horizontal, nbanks = 2)) |
| 104 | +plot_trials(circle_union_suite; legend_position = (2, 1), legend_kws = (; orientation = :horizontal, nbanks = 2)) |
186 | 105 |
|
187 | 106 | usa_poly_suite = BenchmarkGroup()
|
188 | 107 | usa_difference_suite = usa_poly_suite["difference"] = BenchmarkGroup(["title:USA difference", "subtitle:Tested on CONUS"])
|
|
0 commit comments