From d4dce05ad626ebbdbe64f8e28604dc91d510a7c3 Mon Sep 17 00:00:00 2001 From: Sara <sarasdj@stud.ntnu.no> Date: Mon, 11 Mar 2024 16:43:22 +0100 Subject: [PATCH] upadte: approach 3: split into horizontal, then vertical --- .../__pycache__/get_relation.cpython-311.pyc | Bin 5153 -> 5531 bytes server/map/get_relation.py | 78 ++++++++++++++++-- 2 files changed, 70 insertions(+), 8 deletions(-) diff --git a/server/map/__pycache__/get_relation.cpython-311.pyc b/server/map/__pycache__/get_relation.cpython-311.pyc index b6ace0d36253b3b00b4ae6cab83ec4b462051706..bb596e391fabaf33934326971cb68302724bb81d 100644 GIT binary patch delta 2888 zcmah~-EZ6073a4ok@~bG+i_(pP8>&eo6O#tW?h{oaU7>f<76m;4vITL5ouYrzIsV# zQ9w;9Rt!R5piUVe?v?>6kYxs%2M^tYi=o^8fIUb7v1lBOdcgXyrwJGe49m;Tr9R@= z1$IIH?z!iSm*<{y?x#n8Iqv$6)9D~!JQ(?xe8YR!)kogCJ33U#a8ZsVBtjz7whSM& zTQr?^WQ3@Y5u;+p8FjYhu8ceCZquHqS7Kz}0wu8@GSME1`;dtGB^zKs;sJva`3@1W z{}X?4F2b0$tGRR`nai5&)m$p8nvS<qS@|6mc~jY>DPCSltEsle-9KESPBWhiBt@2z zM3L2)bHsb)BBu)2S8HgTWPj{<L=nVi@U}+Lrr9K+zq2HYJgm80YrDpSr)bfkI=f1Q z#%X)~2%6CL^Ic1?C~9^ve6{F&!BEtkdp6+G*s4wcJLTsz7xtik%Jf|Ya`z_0)f3=; zH-}rjzWRYmKi|?4-=Q`3vG|CDpM<wHif;Y4RD|^uy_#45it6ogz9=UN;M?ifNSA(w z9wOa(obDC8l3VgfUdi{kNB<f98ttz7@GkW$jEgD)N%xa}{TUsh15!{QWv<h{Y9A;B z$X;+V$Vmj8WZ&n1WI`nA)3@1cXNRYe99QAw!by~p!ttyWPOPXFBi}k^CHad_@c81* z2x;0AaylJTQnzJ<Pb)$qQO|M%Hj~Stcs8k5xcdi9A%T6$G3?RwfxB4SRV<Z_sXxi- z-?x2X?=cxgPA?#b{>*lyBuJ^7DM^-MnRp&~!L#YW*0D6CmSWl$(d`a-6RB;I?{G~= z=NqvJQ|x>zA4jUmchq<3Oq#5kN`qT^A(ccptIV?wixXQ?a#^JZIN|*R`LL##l+_rL z(-5;<7F__A9bV-PO8<oapmeOx95<NbHRgE3=PzGsa)ke-CJ|&0TTrom9SB!@H>YcX zll8z!BXF|f*zpgpzkly?#op);VAMa{BwThd0;pX0x_^9QcJoK|{!>Q(smjHjfstzb z!BpjKPz)Vu60GwT3)cAU_d>@v#jVA9Xx0eLt~nZ>fv29~C!XQz>l^WhQ+3aT;hCs; zCcrQ}wsG)5aZOwkcY@)n_g<zlTbXT!;CtVM38%*pU#c$G#L2oiX^4~MnQu74?b9pv zOWwt;MD4BX<%@M@-eBfy%)I^&`_CcgmL2`!DCG8iLhWWaL1@HMceO-85e2_jq@WYf z@G=?&rMOgNRU1UGJv8>kB7mMjRvDb)gL_31xojDs)1bm)mEc|#6x8S{hedsRMI^C$ zO=DWZeRd#RR^1;912m_)^+Dl42|5y|*%l}$A8Ic7`!ApV^(*afQxV=|(kVqnZrF** zC2~0=DTu)fUIKnRn~+UDk;|&_R8}#$cs?&@C4?J9gm0xe*k@&v%i~U}m^NiGo|hqu z+XSSlB`INzbO8jdc~ZfRrIVvpW}*?`Snw7Q%I^Tm#J1DBHfs2W4PT`095b9_<?}l( zPkFY10n)(N*CaR=k_e!}ZhHc2@?G(sSYaE2`=hBlQ|p1cFk}crrRvBNVYDWUHrzcQ z74H<+uhrd$4ELcXK|7B$JiYf_AG_99>Yfq9Gg7mjoqf#zKoYqEYwy;bA;TFepWk-+ zpD%U<@7jsVsT$Ld{IK&x_bPw|#8dyN<Jq_iXlOzFjoeCBilah=LLOk53=WNgEyCee zD0=jV!po$j*M&ouvAYP1hG9ZD(GZRt!bOWN0)Pm&8*~N0<T5f!%93fj0S74~1^ocH z3JYL`8u|2twmuHc10$~gP)wEf?G-!6$NTRSRxLpO!=g3j;XjI0kuEZj^1JVd?!Ri3 z#%SatkbTyb!)4i}31aEFF0HZJ$X;a|)D;h}UE`p(#)=NjUZrvVA5ifhw}<A?=miRp zU7(T_z8GAq<C@K~5H$YVEGSI?3kGA$q68MaX16MWoQtp~XB~@TQrSEtI=&2KGzT9< zID8SY73Cn%CYyn4;8trYz!bk_tzasc{H<84icG%Hp*)K@l)8<JDBbS-uA6cRsi{Pl zLEbXSf+B)<>!8}uN`}LD6Buy2Q2q>1CK?BiJ{+twgXNjZt7{|cfo->UwXhzjveiP} zJz=;f$``(c)#&Z7G5$uNw>;M@K{5xM1j(KwTUa3*eDEnBdcuc(wX`v@A#aRrzP~lH zxm1gssSlnt2G7>{X@j4xSr61+%ly1$exAhUeNXwJC;ZT_9UIp+&um=W3jB%Oa@5AB z>)~@o_*|WzG5DF9_3Q|tf%PBPgz*|PZWT<V*R&_)TpsT4xPrp^1K(i^ar&3Oq0&+K zV^gpiOS*6@msj?bf$7b})qFaqrc*bLK~R@vgkHfmY<4A+FQ6EvX=rlRZP#iuxB{%0 zB3zW#4y<7RmwzLZlUCC5DWre^_x{(FzW_FAk|Y~MZ>{q*h?&~|d-OX!SKj5O$cbvA fNx-xhPE$SPfhN&~?{Ik81;3koxk>B?YrFpg5>NI~ delta 2700 zcma(TU2oe|^j_Ps<JxJG#_3j4+Dv0xqd~i=h-qkJpVbCiH_De-V5+j(2`SFUt}~#M zqZF}0k@~RW2`cdtgj5RaLnnAB5)b7E0Lf%QWocqB6PiFG2nL8J&UMl?=`_ZP@A19o z>z;G2&$;<<`$v84pZI)j2;fKQSM9m>nf8r1JF~6#WsFp$Vx7+jaj!dPburTxZ_7w= zDI>?_2H%%ykGD7G{&+y;v|y4^Jr}rmhss|-@sKJ2R#f~fig|C*eNM#KnT76|?~gHy z?B|;?gP$o#lwJo_RNGR_97gC0j5-uODOBW5?l3yNttb?Ht8|1+!Su{AAPb`*>oxhZ zU<ErA#q_qC!DzbA8elhs_ywtzQxvUJY)oh?N~TnnTxC^m6~V}qtgGxU?6dZD^l+wF z4p2GSl+<?k`PINe(Av)JV#RU?WmprO-y@r-9HJ{m?XW)PVr)oNtnav&Sh>8hP7YIR zyow+@ou8*X4`Lj&F7nUs2|cZmWR7Ig*;G_bzZM-!=z5Gv{W7xDbMfbg_W$&I@3q|M z=!yM6OS-UiKM6wyKE*_`x`l|fPk5gp9&0+V(fU*PQuH~Tq3KB?T1(!oC9$9xg;717 z)f_piy*`@H=@Y43)*t~e;Ydv+$5Tk_n&Gg?bc$@G+%}gmI&OeGJ&9;DkRTP2Q(9q^ zXnLZM&SlAdAbDXJ-@<`TqU)S!&57qry;Ux1bI}SHy&ia=(%oO#aj+6OR1F-q1BWM{ zyb%gd9$yJXYKZZ~Tr?$I?~KiFTijjk+-rC4osw3($`9ViHE(40Sk>ESd;2Qxb5r)) za(7v;$b(gR(3S@$hi<k9r}cNU)7dwVOdgu*SmF40hAp4;crdwisB+}g<Wp5HZgcSp z7q?!N-o_#8N2$yDuKo6Q8ea@M-ZQy#Sv5f>Vhr&E$>C_i4EoJQqi)Eg!}`_#Fn-MH z4s1RhqFNrZ30{P!NT7(Oz~S_?Q2?HsXIQ`i@EMJyG}RH#<cO*fgX{)D!vn}oY~B23 zYJs9b>$5<*M0Xgw^U1@477z>LcNy>p$b`!6f@K}XS0|{*0Ebc0MP`i)MRlRb7N~Xu zOr}}eWY_c?^pY#sT6gOVG;PpgpozXtpGMpOaDZKNutndc*Y`*NX@I+jfhKyZev&IP z(-)!9a3+GBd;>Y5=hlF6ZgmFEm71b}Ws^x0&p{%u+({-FnlF0pX{hOt*2R3cowuL^ z@}>s{E}aKwV2bG*{1%!Us429Zo}e~4n`D67oZ9<wn?iGS>ke;sa{;DW*eR~HV-~$8 zTjpt->s;T7N!aA-O}ez61x@<G^rzF`7HB7Fz__?NGr0xNx?Mz5gyzYgAhaHc0_YZ< z-J)keIswJn9NsNp=wi@M0VUpWm~7tR6C_c}rZhrDowjtg0OhBR6%63y^IBFVz0}4Q z*M?h`jyR_2`lyjUuQ~p)T;>cEusT{d?$$KfOHJXlgU21om2K=SyNgDjb;NoPO{vL} zLGNC-($~t}@k7f`tJ+!*jf2o07{-@w553XZGkf;(t9EDK<Vk3fvj=7m&K|tHr-nGs z5nRX6E+aj&rTMdV@6N?tUp%qA$sRmhjSSh5p=L(uo}B;3s0nQ32%hp!`D-%r`6q{J zJW{%A$nP0&(UfpA7`ogwyKjD{9eHHFv~<=UIIz6C8XUHR!!_h9i6G|_2c`~O4{x}9 zYSrO=%RBACq2-Bcc*G8m)R620F+i8_joxT^Q+Z^*>!Z_)2bRRE!<C+cm!#V~if)^a zERHPnEJc<_zUry;9bJh%djEB)NdaQ%dkX)C*yu6*hl?-C5S0@7A++UMXiHhIh6e4> zU=88mz>3oK!SVNwm%FP<zpeCFl>WbKUDQh!g1H?<LB$S4=eJY>JF0;lc3=nGXp&Ol z6yjbEVMk19xqKq4atVX<T7QKevZ6|NX)FDDB;%DtUehP~b9rO^{SnM03VA(O(9>u7 z3+b^}H9}}HI-cygOn!o#q;u9N5ckq`FAPFI9YQZMcWd`L`Nz8_lT**>+7o0981#RZ z@iidu9%H<MHdI!h6|}2z@5j2WoIJ|);hp8N8iLn)^aAGb12uF9b!@`n8oC>~p8X4o CXnF_$ diff --git a/server/map/get_relation.py b/server/map/get_relation.py index 25688489..9b179982 100644 --- a/server/map/get_relation.py +++ b/server/map/get_relation.py @@ -1,5 +1,5 @@ import geopandas as gpd -from shapely.geometry import Polygon, Point, LineString +from shapely.geometry import Polygon, Point, LineString, MultiPolygon import matplotlib.pyplot as plt from shapely.ops import linemerge, unary_union, polygonize import matplotlib.ticker as ticker @@ -21,20 +21,41 @@ def get_relation(self, body_of_water: str): print("Failed to convert to polygons") return - print("Performing div call") - test_line = LineString([(10.4600, 60.7451), (11.2000, 60.7451)]) + print("Creating grid and cutting polygons") + divided_map = [] - new_polygons = [] + # Divide all polygons into sections for polygon in polygons: - new_polygon = cut_polygon_by_line(polygon, test_line) - new_polygons.extend(new_polygon) + # Initialise grid + grid_lines = create_grid(polygon, cell_size=0.1) - tiles = gpd.GeoDataFrame(geometry=new_polygons) + hrz_lines = grid_lines[0] # Horizontal lines + vrt_lines = grid_lines[1] # Vertical lines + # Cut polygon into horizontal sections + for hrz_line in hrz_lines: + # Splits shape into upper and lower section by hrz_line + cut_poly_1 = cut_polygon_in_two(polygon, hrz_line) + + polygon_part = cut_poly_1[0] # Save upper horizontal section + + # Cut each horizontal section into vertical sections + for vrt_line in vrt_lines: + cut_poly_2 = cut_polygon_in_two(polygon_part, vrt_line) + divided_map.extend(cut_poly_2[0]) # Append part to final list of shapes + + # Set polygon_part to the remaining, un-split, horizontal section for next iteration + polygon_part = cut_poly_2[1] + + polygon = cut_poly_1[1] # Set polygon to the remaining, un-split shape for next iteration + + tiles = gpd.GeoDataFrame(geometry=divided_map) + + # NB: test plots fig, ax = plt.subplots() ax.set_aspect(1.5) ax.xaxis.set_major_locator(ticker.MultipleLocator(0.2)) - tiles.plot(ax=ax, color='blue', edgecolor='blue', alpha=0.5, label='Original Polygon') + tiles.plot(ax=ax, color='blue', edgecolor='orange', alpha=0.5, label='Original Polygon') plt.show() # Convert GeoDataFrame to GeoJSON @@ -49,7 +70,30 @@ def get_relation(self, body_of_water: str): self.wfile.write(tiles_json.encode('utf-8')) +def cut_polygon_in_two(polygon, line): + points = list(polygon.exterior.coords) # Extract polygon coordinates + + shape_1 = [] + shape_2 = [] + + # Split points into separate shapes based on their location relative to line + for point in points: + point = Point(point) # Convert coordinates to Shapely Point object + if line.distance(point) < 1e-10: # Floating point error tolerance + if line.contains(point): # Determine if point goes into shape 1 or shape 2 + shape_1.append(point) + else: + shape_2.append(point) + + # Convert the lists of points to Polygon objects + poly_1 = Polygon(shape_1) + poly_2 = Polygon(shape_2) + + return poly_1, poly_2 + + # From https://stackoverflow.com/questions/39338550/cut-a-polygon-with-two-lines-in-shapely +# Splits a polygon into two parts based on a dividing line def cut_polygon_by_line(polygon, line): if polygon.geom_type == 'Polygon': polygon = [polygon] @@ -58,3 +102,21 @@ def cut_polygon_by_line(polygon, line): polygons = polygonize(borders) return list(polygons) + +# Create a grid of lines spaced by a cell size +def create_grid(polygon, cell_size): + min_x, min_y, max_x, max_y = polygon.bounds + x_coords = np.arange(min_x, max_x, cell_size) + y_coords = np.arange(min_y, max_y, cell_size) + horizontal_lines = [] + vertical_lines = [] + + # Vertical lines from top to bottom + for x in x_coords: + vertical_lines.append(LineString([(x, min_y), (x, max_y)])) + + # Horizontal lines from top to bottom + for y in y_coords: + horizontal_lines.append(LineString([(min_x, y), (max_x, y)])) + + return horizontal_lines, vertical_lines -- GitLab