From 39f150bdb7421460a11bfecdb57ccb909396439f Mon Sep 17 00:00:00 2001 From: Markus Metz Date: Mon, 4 May 2026 18:21:21 +0200 Subject: [PATCH 1/7] robust signed polygon size --- lib/gis/area_poly2.c | 17 ++++++++++++++--- lib/vector/diglib/poly.c | 16 +++++++++++++--- 2 files changed, 27 insertions(+), 6 deletions(-) diff --git a/lib/gis/area_poly2.c b/lib/gis/area_poly2.c index 3fc59eefe8f..1cdecd92e7a 100644 --- a/lib/gis/area_poly2.c +++ b/lib/gis/area_poly2.c @@ -24,11 +24,19 @@ */ double G_planimetric_polygon_area(const double *x, const double *y, int n) { - double x1, y1, x2, y2; + double x1, y1, x2, y2, xshift, yshift; double area; - x2 = x[n - 1]; - y2 = y[n - 1]; + /* get a reference point + * do not overdo it, the first point would be good enough + * particularly for small polygons + * shift coordinates towards this reference point to make + * calculation of the signed area more robust */ + xshift = (x[0] + x[n / 2]) / 2.; + yshift = (y[0] + y[n / 2]) / 2.; + + x2 = x[n - 1] - xshift; + y2 = y[n - 1] - yshift; area = 0; while (--n >= 0) { @@ -38,6 +46,9 @@ double G_planimetric_polygon_area(const double *x, const double *y, int n) x2 = *x++; y2 = *y++; + x2 -= xshift; + y2 -= yshift; + area += (y2 + y1) * (x2 - x1); } diff --git a/lib/vector/diglib/poly.c b/lib/vector/diglib/poly.c index 74e1622cc70..c13ed8c84c5 100644 --- a/lib/vector/diglib/poly.c +++ b/lib/vector/diglib/poly.c @@ -96,19 +96,29 @@ int dig_get_poly_points(int n_lines, struct line_pnts **LPoints, */ int dig_find_area_poly(struct line_pnts *Points, double *totalarea) { - int i; - double *x, *y; + int i, mid; + double *x, *y, xshift, yshift; double tot_area; x = Points->x; y = Points->y; + /* get a reference point + * do not overdo it, the first point would be good enough + * particularly for small polygons + * shift coordinates towards this reference point to make + * calculation of the signed area more robust */ + mid = Points->n_points / 2; + xshift = (x[0] + x[mid]) / 2.; + yshift = (y[0] + y[mid]) / 2.; + /* line integral: *Points do not need to be pruned */ /* surveyor's formula is more common, but more prone to * fp precision limit errors, and *Points would need to be pruned */ tot_area = 0.0; for (i = 1; i < Points->n_points; i++) { - tot_area += (x[i] - x[i - 1]) * (y[i] + y[i - 1]); + tot_area += ((x[i] - xshift) - (x[i - 1] - xshift)) * + ((y[i] - yshift) + (y[i - 1] - yshift)); } *totalarea = 0.5 * tot_area; From a00492d24d0fb9fd59d06bd03d29704108b7fb7f Mon Sep 17 00:00:00 2001 From: Markus Metz Date: Tue, 5 May 2026 10:45:19 +0200 Subject: [PATCH 2/7] add test --- .../v.build/testsuite/data/tiny_polygons.pack | Bin 0 -> 15875 bytes vector/v.build/testsuite/test_v_build.py | 37 ++++++++++++++++++ 2 files changed, 37 insertions(+) create mode 100644 vector/v.build/testsuite/data/tiny_polygons.pack diff --git a/vector/v.build/testsuite/data/tiny_polygons.pack b/vector/v.build/testsuite/data/tiny_polygons.pack new file mode 100644 index 0000000000000000000000000000000000000000..167516a8c7be1ab62079854eb60085df40518ce3 GIT binary patch literal 15875 zcmZvCRahKdux)U6cXxLiT!Xs=clY3f3>w@aI0S+_!QEYg1c$*waCf`&pZjti&P#t^ z_0%p|tE#$Z_oj$NM7&#jw1@uqYVYD_>gMX?XY1GK+(;l=oD?Jh>T=3Z%FGZ%-&ktj!b$2;*g8Btb*QJ^ zK+4AI^76vMu?zpryVBC1gSE9gZuyb56_v|%+b2~W9UaWV^ui-zLtH=W7LpHHd0v_H zSbMeYovQg}P++!59k~oqOPL}a+L$Lc6S(hv2TgApX|rNBBZgj&}u3H}xQ`jGfRK=^@!js3d!Ve0Za%QlsBeR%s3 z-X^&8NMRcl|(?}bkeQZb@F?}^i89usz9r#Ylk2K1-xr6>tvXDWw}a$P+17-i!fsgJdKDc86uBz4V1vqjbXp z#UdF=nGt}LkO0Yj*0Agh^*Mcbu=z6JaoBt*a0P6>9QX(}UlN=J+mqGTga||d+2Z)3 zgIwt@ond?O`aF>UiG4npzG9#Y+Dm9?4%$m&SQ8oWFli_3x0V1sh|)G&Q% zKn1jy^iT+~+lb4wZk+!y#GFq?yN{U0_9;jX91o3-e<=pSz;&yCG$!rS3fFk|Z`NH9XBH?vEstHlocPCVCli?BcEH-za^KlvRqcNiVlfE7o2J_7BC8 zWt0bb)BBt+x}o4i`@U123B;DA)39zI-Za|V>!&Eh@6%_kM%lY(eleViYdf@hna0l$ zyO&z5l2Bf|KJXM?Xw8@>`yKaqnh=Czq3?h}9Ekhv0(Jc=dIZ)XHD5gH{AA{vxu@+} z)a^{1;&`46dP|7d`bTBQ|#ft4oZGH@^w|+e@Lz!b#f3746cd( zb_SF@UdDL3!VMx2{d#SXP}p`5I|Bd!?Xpj^Nn*awI}}u#m}|O@G@BklulVz-vPh=V z&&Mr$ACbJb>3>G4yjb(UqG|$(-b+Pl)LJYt!fcSYtA`7@y20knC$$`NRjARR?no01 zXOw8Z+V!B=neb@@Lu7smCnjQ?|E?>*`k!{nD?n8lI?>UJ8$4BZ9;Us`6;Ho@)f>2Ynel=pqA8QyQp92cM5mq~ckj+N&= z&2U<*1uD@{TK8gPa$DTwKLl&XC{G&JWdSpA`2J3xoCNDnAQTMtn)|=Gf&zl+ly^jt zd;VkvnS9)tz*M!UkhKI6e~Rbj{)>odni%7&WEpxloi_2k5dR0Az@As`cG}Od)}SV& z`3Ev01|>qCmju)FON(64SBIIscw_x|R~ND7hUm6-uFreA_e49pmb)rRwZ~g7DJZcF zJlo!-{khwm<}>cohh1yo`-Dwq<9EUfgT~^*DwDZcTjDdalXkQm*&_E~lA5vD8iTFB zW!P5CO@_U$)8Dj?>&OsjhHu_n#o&g}MQRX~E}bbz(%s+8x4v9&TOO586}{~|#C`<* zd-Z*f&Q>uyx(!J1$qf+_-@Zd%fEo!GC_DZ5`RV%=J=gdRBH`3`i(M8l%8{!pZNb#@YR*NjIQUIl zpv&Td3!Oyh_9w@<*{66ixdb$M5j=_=zKUsV)*1te`82ExU|nL28UUPs*)XAo<;Oz zd|kE)dV;T9_YEAwGB8K=qvX#<8~b=!M)B7e67KciUA^_)$uJvyUmd@3TQ<~J>!Zda zO}1QtNk)6knKBO6IyRBLN}Z>r=R@>gQq@}W&Frlv@l+{Q$^Wp*+h)55Da6(ha<%wf zczsyd%M2iVK^JX2gFFJ%T{qtw-ZMcVcDb~+MeZB&9a*YzNcI{%AqQ5@e+SO7Y|%FZ4pqaqyjbQh^U zig`NxawrcjBjA7CyWmt7STQvhpP3g@k=sBNnhJc&ZuVDIcn$JVz%Nm#e)WA4&O?fB zt?o33FHv~)*=s}iN>5uWJa2ME>7V?2qP>%Da65VG2$Selw&EFLLG!6HX~Rx(D#YV2E}hzDJhI25#qLHC)}*o>dm{3IBfNrC0m6Q&-f$W@2)WToRQl?>Ok=-z{ld{LVGn z@yq++NA*Rua9@7+%BNHu?MNDd+BA~ae%{zsDrUGpClZgp z4({l7!^-(!&3sH%#iq}~LCvJea&T0=`kyJcs^k|(E&VLryo&~EIosa+CXiH?b**0opH}=qNg1$r{x+TUDFQUwg;Jd1@92EPyDtIPF zp|K6;oFnDp;4r#yr=<5Qu_1J`y7({IMwrnGkcp=TYhjw_6Kr|6&ZG%{tj$ms>o}2L zAVV#Um#&aWoVEx2rp`MW+_M$|wXreVGARVFV-sv64Z(??4GS+^U0fS35M*b0eEl^W zOE|t>c5Avt724(1Dl)OX&^d|&n#_--+-|b6EEW~p{;ePuq}uK7Nt~Q$^4ohd@Od;S z__>1(Az}!_em$}ASW5W1v-B?fo)LOJ45NEfz3wG9-B`&4Kb5D1r1-%<-0i7SFBZ2= z1SSLKF}!{&UMzwPIx43~xSKlJ^ z5;n9^HMs#JnqMA*%*#mY`ULY&t~y`Caht!coq%EaE}A;vG^}(7qbrwECo`X>TtHUH z0HHlq75nHLdew-t2ifWytOn?DnHICl$fN|B`+g+ne|b&3&5i6ZaM+SZ!}*pDK9s5O zGX`pb{9I$UDTuTcx z&kR3F5=@+vLLt&fg%PEnQcI zgCnSp2QPkzM_xHimT;N&!Vot58M12?0*!=EDj22gaARj$ z0flj@Sd0{Qj*$s8iE(Oc-KZk|lkJt;ItC(e|C5K8NGx{3nGK8tixhQ9RMaDahImsbS3zBYVm-& zB#d;s%rHM?deHHI(A1i$e;v9*(W65=7@BMFN^F6P&n28v|IBsK*CIEk?T)#ZjYf56 z&ZT}lD2tKl6RsbDjT@uGC~2?Es9OXZVS%Tp4!mL=Voobmua-hmw3LcQil9mRt}eLy zBNj7KoiBC(v(k=_K21*jQon=XtBDN9^@zzIKM${O^HelYzDoPZU7m;TuAKc@1N*HCg z-(ydENUPEM^Jly32Y7VFsIPiam@p+&t4skBD1gUWI73RZK-OI*g|MFWG&K9Yd&}U6 ze)}mvCyplWaNGQse?{(OAb&d}c{?1I;{uk*DL*v1 zu`HKVEbFu4NKMLN=^b>dI2dO!=0 zuSP&8YM!Ahl90E9)81>37{$uBR|1}g`Lo_#Hh10W_sMK<>qlk6Bi31Wpwq*?TdFJ~1Fi#{SmaIBWSPT*DR5Ho9KeG16>#R*c?#CYX-qVw-Qf}~jU4r9fz1UIM z^JbyS2~+8u$=mboB=?=-?1TbW)+NK#{bkeCfk#y1XH)s_%3&-UKk;0&g0iYsc$wDt z#tn4$zgJ$Fk<)AaX+zfXx6c|UcH|x{`;MbI^_gC&noiPk@-#~<-PieAW}ZPXdoIarxfz<=0EoGy3G%DayYUyRGcky^P?ui-M1jx#t4o48h)y z%wF)Tss4E_1k0h9xAg!M7&Yxae&~wCDR>R%e zx)-gVu@WCnkd>KxGG8~=LuEgM%$S(1oF!$g&_(pgh^7>~{`BH_otb@@EEp2;w*5AI|gsiKayh9>3#D;r1s z5Q705eeK$3uJ*hGm-dyJYma-^29{ZaG*Kr0ab5tArG5-ODRPk-)1~ZNU7h+o@W$f@ zi{cHbOIgMDz}VwzG_BgI@*i!qd<=#G4&I+--4mLocpC-``1uUyNPB*%rj80t=5n{x zZ)gwnSbUb%&?uX1qE4+UqbBq7d+pQ-HogwB=J&$JJJgbVgSWQ69VfhBuFdV@);s!H zn}>4WHrc4l+OVX>Vm*T`-*BV3GzoW7+}Z2Fo*7*+mb}NBVoRH_SVsMR(^wEzSvA96 zEc9(zS$AZtdh(azRB(Tmd{tWSS1%)atpr`H&9^0YrT~Os7a>+V((`knwceiCob!aq z_EG(+q8NRwtgXUwx9kP&yf0|R*5D-@H)}Z!^XR64gxElL`PId9sphLPW}XdeRa?K3 z=!&(G_G?X=8zVfpQ=w|W-bC)fo{_9yMZ<4=t>j5Qs&y#Cq(0_|I~hmqGE2liOW4Q%c=l%x!UFtohuIJ>6bZ`*+V0|WyLbu z2>;?_y->fVe}``b&EeX{sHjFODKw0H7%aKd?JHy}w>KeXqd260zB86MszQ?zN^71n z=<~P4o6zU7IC`?VL3l2=Z>Flh#r_57=@COGHk_?JZ4B9LiGCJ%D)1gd`P!>=9Z_@o zF?5{T@XLLSHrP;0`+e#CdMo)us`#y7VPs#@G#@o|@J>GExp{azGEZ#;aMIozeG zXCB?zSy!KV!s?#cj(lU5v>49vPNxxb*xS~kAlj6DZEQV^&VQSM-_ZYd{Qhi2@pUIZ zxViUn7*5PV`jX);cw7GUA+L8s!)QVmRUKawg|2<|;OFq_lk#eo=s0xjC!PvKqYaZ1 z=3Q2%D~pq1*N=nNkL%IkTT{KtH@zpiZy)J|1g;60d%*^lY(Xy^(q2y==<}hAgkh9m z_@NU-tSf#__9cI;SQ??|(dsqYAcQpkJ#Wy%IUXJW4KqtE_w4<+(F*Zk+@mYagnLGeJocOxRu#rfl`-x(F0RS#U-W^=38KKf$JM*3iNi?BeML3F2ri_Yh#{TuWY8pKgB6G zEMUOyH}cH{DwT++mDGRP*v1;iH#N}P$-b2`80ir03(E7%xif4?yiY@`bRBk5ki-tH zOGSeqCG<>^P{(nZ}%yh;@N zoTnVR=87B$btonInTQKtqG$Tn_f0#*{?A*4vn)727j^{Ii z?M)NRmn8`EH`N2{VT%9y7p$9^h5Kkd)J)aED+7T$CtSD>6T#n5100a_50#Jy%AOf0 z9NihxLkQu=7PHRh-x)fDRM84_OxJmFYoM`hVV@ek>f`SxF~|}oE41(iPEPp99qPv&+64HKS&3}q0@>Xt zGYNEvlPubD1pfqCNU6;CzFrw@|6nk=?Q?8H3?f8=50f;z^{R#z4Y=0Jo(O z%4s$S_neEeLy`)vAp9;#A%gqphMfa#sX?g{;VL1eWrm0h2@#I=9VfDE?ob#m=#{SA zFfQ2vvq1xr#WF3G7MFGMx||AeRF9q!3aP-bgRY zGhb?cRz!f=7>WhAMLz#*iqqEbS>5oV~@`M?EHMt>=nXBoPS zfENFUagJwpN)*Ru*@oDI6>SlBWq13e3UXB#FtRMM3{=$uJCbx!xp0P~K`0ld_2q_F zfy@unWpI3LJg6Lkd6Y4fW8o)m@i9;28T;!Eq2_jfpjLePprQYI=>|f*_;}#l#@dUs zj@L_$tIKs=uwC!4`Zyu=Dr62YU{fN(z#Ij1A_IA&TjC_czlvr;Q%F-=n9-tkJ0SyE zk;p=e0JjZ#AP%bEtX#-(wJ^Hi*t(>SVc%}TvSo;c;4pzW4kW~CV`iZEeH2wN*CGcb z%pO{7PPtEi)?=WzqhSB`O8}mU6+jf9EbG1rMkWv6av6Z?d?HS_E!3$|rL5Wz-|?}5 zsA2I(+>y}Y3*dNF(g8KferYbFY!i669>j+`DqOYAJGjQTjfl%T`&SLDQZxw6 zU&91$8Q}FH{VM)~{g0GOic$`i7eY^%4AD?|sQ-sHxc?oS2HBF<98wrmv~!%1u8@8( z{a4@+IARqi5t%%SO3**c?F+MQ!N>zq5g_T>nd3wyW}A$-{23m}HMt*P)inu@_m9=6 zghVt;)iBQ~mo0YBa@lGE^}PRynrECFegdn}=vGnW!iaUYs`BXwgucaP!goGyc`&8A z)gZm(dcb9^P2iO?bWgeL3~a4x-L}vciyN_pG$IHzB!Qp&0N2}VpEFC%G&IZGflR|o30$6~0)m-*rynIMO+$4n8 zTsC03tWro)L9MMgImUbR(ko((c(#f=OmN?TEZv5Gp`@NVsysiV*MJP)^f^x2N?>AK zl>O!ObQymrTfoMI7!9LJ#050J7PYS^d<~>gf<>aQVj7I^+$BJ!Q}3MQtRo-@SHZ3Q zB=q=$gOEvUdX3_{QjJ zn7&Tdg3w-*@VV{zm|tHws$5uru|Um9WrnJ9Msng1#}F&~xa()Go|2SSo?6;XB_(gV zXkG3y?saYDO~!SrfyA}eD)HoO@R0f&f3;m|O=G8xHP8w012`3Hd9h05>|>ICBG7-IhzTsFe295`h|-(1AN!bR zcB6}dPjA`-FN_5&64}SV z)Lgv+=7}2S8`1rj@)_=lZBM!ige4)R_ATAyZoefW@pTghF-DQd+Bc<9?e4vh7z-H|S1|+Lz)hqlhUJ9|aliq!1?s`(3%c zd-+d`8H&~%bO-4N7mZ>8H;;MR_n6l{aKXyUWJHE2b~JzEo`!ANb-|8FuiD-Q*EgPM zAKja>rMbsTt!` zHsYn+Ff+7hGv^SmvRAzYn@!TKTYFJHRR9wwW-k<+Dzxu48@p$vGtHDEeQGCyv|0kQ zO(I0gC+HpM$JMNFn8ZK+GDAj$oGMO0EkR?m_NuiC-F+7Gx`I$0kxrVVwc3*%o@eyL z`CY?@StXPwDj}%foc1 z0+IHi{6Qld^&6T)*Yo%vJWmDN;l^{!rFh-i%M4hlFO|R=t zuv?^p7{=E?JZtp>FH^jWO8?pBl8wqH#p{63v!|LOB`X>Lv!8gcZs~h zf=GE{xU{BfQq?k{S;ylxO<4z;LDieAW6pZW!`c$o9vAe;6X#JBEvxvu-sy6iP*>TY zNp6ZVTK-aatZX3H@{m)xR%wLbOcv~sp%R{0c?@60M6Y%?Gtj=`)q08cK=-J3|A7<|c7e8s$>2Nq zd-b(-Xi(F>$HukD`8X|ax;T|k9aVVgwIS_c`kTmkmb@LaxUz?d08uS!CyLIx_S26F zQ#h3aN0v6^?1^RFFlIaEUYq91!Z-3>8|GdSDxLkZi7KjbwR=$O#gvk`yfHYRUMEhV zvcU$QDE0Sj;snJ)1dxSF)ZzCP@7W94<1ti$ zC={Mp$RnBd6Mou>1-hSJdLrA?`iA1Q)>A`Df$0OrF{%58#*=UADAVXgn=N}$#cT7p z9~wQcrCXo6qDYtZ%Wsq?^v)V&Qy1o5H#OqgnI;n`nfs&z>?Qp`!C8^*%D)VLJ%M7| zA}z0Fd%siMw|YvCTRlvM_m;e>95CKuQ1yPgI?*&tr=My3R@H8KQl7p%f2zE9J$OsI zHQ(Z!i}<6G$5H=BICIb0r_w)VT{;T;qvw^PArtwHtVNBgchKeMp9_5fxaX=RTeREl zHCm$JQR4EE3;mqMSF#&_PNFF7n;L)^EY`ZhQWtw722`fsCb2A0!*=C!i&81-_A~ z_x(J7%iqw8)0@ND!Cmu@UPE2)HLeMImNd{I>mF3Z$NGanxzV#`rHZZu)T}@p0zPip zGGEo}?eGY{pXVpv8BTKbfATi1RU3b9fZ1I{p8%#j2+75=m0x{F(5o#x_c{Kou^oPRZmz*m)3aSZ`+I$}PLA9SXoeh~boh>L#yuUcWxSACF(WBJq9v(PGQ~ zG3$H37ql8 zN6b%-w>q39uY(=N*H9?Mg12BBpoXD*z&%Dn>tz25$X2L0nAUedyvj8*td1nYoC>u!>SdkT z5Y;}5o&*(oHi|WIFlqT11MuPVa>Kuh`I|A$p+Dr|v`@gN`5&QWW*#;4EW#Z8XOXn{ z!1#eepI=}hPyRV}l6T^WD{w%Zt>*w4xi3C(f1!Jv2BCOHEcf=kbdoWD3T(a)QDCy~ z1KiAWr&X~+F{cZwC9vxJ*>^dYhrg3_V}KDA{v+`!x0>hoLc3R*F}e*D;DS?Fu-EtR z@1zNN@_hT$JfBaqR3k76&icJqyiZp_Ij6V?zVriKpKgAK1&Agko{|nRqCWOh2LB7U zBtP^_T(`0e0t<~@ua3Uk-2OB0v)XIXUiIB(M7`|?llxgo2Z7xi0O<9-_J+GY{khZ# z2?8;L)NYT1PHNR%E`NXUKB*&Jb`4*zVBE}fjpP&4ot-uon-mMKa3PN*BygqI;aAs> z!>jh)@D~qK#qz*xD+esBEwEFNL*;9Zz-U_UoSup77&L~Ng8fpj`ub5AFQjLGt{uQ| z%emK^aqN+Nt=FOh+CF`PzI3wv>ZW{u8Gzyc+zlL|NBG8zd}!v%*;Dazlch;drFC3& zD8v?{c+}B(Z^w99CAw)6)DS=jlly`*duIkuW)#{g+D0#A_gHp(?IOQESLk=)X;#VF*qBzvc-sgJ}3#?PsX@{nXO_?)LSu??eWxm zYQgO!<%RQ&_=vnp=m?$&7M!<7^eB)%4Ip`2>1O~xFT0SixyL7Q(5)=cnW!Qr{ zq!7FqgGvI~lI$JH`tJ8l$nBI{AxuqAcWQ6l4n6U%79m{%X8#fW2ws|fgg~zamyZ;O zfDT&+yxufcw-Hir)DR?p64_h#w9UO|U6qX15Huv3^@R>@(b-yh+TmSoN4j*f{_3W{ zZEyDW>#i94v{4!|Bdh@6(}QAL-P12^LI`VL<8f#6RHRQm`h3AeqrQE= zl@l>+?Jw97DE zapeMX2)Ef-1HYPvGL)P1bozZo=_C#3Q71f6iMT@<)8|9*3!(5m_!I^j7{iLAxd>~v z=?==rST$;Yv4Ut?!TGwH1I*mo3@#JiNB#;8JsuS1g9hA(h({1q8B#pdPKDQ1&*?iPHyD^wzAxf{OK*U9AA-w!+$~ z-LADNVvlCH2#BOo`JF<4lT^g)AQF7O8Ne_fP~FV-WCnF9D9LYT|EhH1xANUQHILT8 zgK!T@KFSGDvGIj_gI?tG;>471UxWFgCsz7>LYPydmJHZ0?ka?4@Z>LlCE&)~HRlu- z%a`050Dg489iva5t*&3Siu%@fgSZFPt?+r-CVg<(t9~54BjnH=@*sH35bgZ!xgbrk zvLABacZ0bHW$U#(>$}x-e;JfR+9#$?_%)R3A57fClXh8I8~28ouLeygDJ29@LlsNl zLsRc!3=#Q$W|%iTaufkr--68a!7txw{!fsPgM24x5|g@M$M@mKKaIVwOz!{bk@nl1 zS&i)BCRpnBKakJOE&OdZo1N_a)g(=l0i#beT$BmWzqu1R=gRolq^D5kiV35A;Lh(V zwiPwT5Vc9g$c05l5C9W0L|yTTn9eR-*4*C*zNDOa6xK-W@Y-O~LI} z{e%=7#VS&-25BW3&GGK^MV$BQQI7_xOs7&B!j_VEBL={STUgNIbgKzewpD%a{dMeH z0s+XDY7eR`X<3hT6R%UGKOuqqcM;9y1{8PgTwR;6Ky&e+<5wp%fNUHBso;#unmqvk zg$(P*o)|c5)|VV>)_-AhJRt$~cR$tT#+a{J6nkZP={@UL|M;1j_+83Fg^`EV?Ew~X z!OJ2Qu@DPf1c5`)hOa|wS;>hZzYu4X1}Swn=*nuCh(#bLzhxXjG@cNgaP|6j+L`Zp z*1EbMRw?JMey=+~Ly^rFBVF6dlaF@x$9Et|2=dBP3d#HTxhCVK(YZ_X5}LYUh*fMB3$cLOD4O6Uh!xcJ+vw1VYT+BK^qg@{IX)WV<=l=hXL@{ z9E;znXW=}$fm+u^Fl>@W?rXW!D4}ii&FvDOsN9OcPKwF1uyy;2K4{1KrVTolg{!XP3a0O*6OmJbZ2B)hx1d%lr8@ zo%8EFjvh1TxI#ba^=Q0G$C~+`!hSGa4`&;V<0vJY{@Ev{MUibm{p-NaBU_5#$#8k~ zD>pT?eV!(nS^(bD`7^g{5cbp(;_@88adFo_!YbzUOYC;OVse<{6~eGqNhdysCb4X7 z?QovL8q;<0mklK=jF>pYo%1?>4J%E=2g3r@WJ{^XnrO~NTd!fp0O``|!r(RcuhfH1%wkaM zrd`S8lHHr1460dUSgWU(S(df-|E_Jobj4F+EdH&X<6*ikTLbj&;I3Ap&${Bu-j@%G z=58w_!+Etio^px3b)!FUhz^wXF7G&ww(A8#Lojjm^IjxNkMD9puX2t%e2 z_pVO`V@N4RB7S39$P%{cQ;urQ91aNnkf1boLj#bG9E?y;JUK$Xkto?EuPq(S>odvT z*G$UJa>uaNPGTM4w*O+6Z)5UsG)^cl9K$&F7^00qMLaaMt`$|FUZ)) zn;YCxVUfL<#4;dcK^*u%Fd#I*M30B|@DW!kqV>u*Xc<^+5)f#{MJx>rrgr-(vU;(fKp6 z8x0tfM><@;HmKz+j^eQnEWop@Qp_(ox7HvfCJ!sv-aGP`T%vw5B#yE%1e&%o@DlD* zFxCdzD{)NqLQ)mipXVugw5OUAM($p!-aD}vqq(H*j(q9@m$Bg=pO~Y7zmt6aT*24z z;;uAfxx(_;1<$$+$$oe9f4P=3TCrD14j^{VVlSXJ6!=?HOD_AD|8JgI_CG{o4;{-b z>k%yR5KsWgO@g`%`U(PZ&gH$B7OsioO|*)j1%Cr2605@iT_OK%7_0q&?00nBr&fBS zjc-539&+5qsUxW+AGxE?xX7Lqw8VBym6GcO4m<;>WrG&dW>-X+8qImcglyHVJb1(~ z>I%mirXRKWQ>a~?M{K{7FnUPyhz-V-j&)BzqVi-=_d1R+q~hbQmIxInQ!_P|3lW$y zJOgL{0{44X4+HoQ2h|#8ADm_wUJV%=O;E(<+Zeff;3$~SOU6jPTozeS$<8?RN5#nE zKDnWYF>k0O2Yk7NDkPB&+6;XFr;A-{Ds8Qw-z8^$xlTymrQvAa*IH$cA^j6GneWNy z;laWlV57C4plejxnLs^V<3~Y!ys&L}W5AM0De^5})AVaG{Et+bZ}~*jaUp-ddFg!D zkKbUX~;j9WQuZ`a3<NfTe(Z`L`%L&-v8cOmq}gCfv!2KWi{wdIe=1|qw(h_< z7J>C&Si}15%Hf3Ct#PRz+2)U=)lJms&L$Z7eM>(Js>`Wi3C2nM!#zv723h!Fw?DX; z^CFH!n4QAOnVSLNW6V|-Rq%7BIW_y_xSZez;l-yo%viDLP|5oQDZdc1g>XR;l7*N< zGx|eLQclpROMx1Zab*fupy3FeGWtqmN(OZ)uE;#~T^ z!NnPy?ZWI4aV9{?QX_sU`~*Kt!^miE0w>GP$SZk;gsiDv5;}l~+r)_##)V0U1_C92 zo~kwb1xw;n-;Y3oR!;>*;#RMLKtfpm8;XRw-W{5RwH|9hn=qXZ2NL`(T+X-pRYW;f zvs36uW;1;_czv^SICx{TN4Sa4>cpWQSZP^>x3JcI#D>_1l&;|-Zuh1(7!ssCBwUm$7sE^mtrGR&ueo|)}-A7$V z$gNRR@Wf;^dx5_8KkxVG*|xDKjGEb#Lvid~gAV~iZ0xr~@yfjkQsEGt9;R&_-BY)m zY9~odOtbnNWvk9{${LLK9GnWM$bUJmA+pL6k2Wwb+rdv99|yg!nIG4t@00&aAn1;S zU38kS&^cDy0p zDIgk{jr_oT_j=@lfBwHszn{c(hB^?&L(Rg&H3Y|oLC4Dtiop}3&BDUf2MZYC-4MiO)kO)oTA_Sg<4!78mJv;WVO)}F2SF;Kne3tCm6WMU0Vb|rV>>haW#6lj zew*nYck_luVD0{>xs6#sY=|?9SB=A9^GNSYEFapRWBNopM9(mkUCfp-S8>V1|7$2c z?{Lai&GrL7L-6LcIAaI7`Qv`}flZZ vd6TIWq}tlViVOC*DhB-0SJG6r`LF~3A8bLWGu~PLfRT#CBK>3s_38fsPxUZY literal 0 HcmV?d00001 diff --git a/vector/v.build/testsuite/test_v_build.py b/vector/v.build/testsuite/test_v_build.py index 03915edba86..3f322ef7cea 100644 --- a/vector/v.build/testsuite/test_v_build.py +++ b/vector/v.build/testsuite/test_v_build.py @@ -6,6 +6,8 @@ class TestVBuild(TestCase): """Test v.build output against expected output stored in vbuild_output""" + tiny_polygons = "test_tiny_polygons" + @classmethod def setUpClass(cls): """Set up a temporary region and create vector map with test features.""" @@ -64,10 +66,25 @@ def setUpClass(cls): 0 | 2 | 2 ------------------------------------------------------------------------------------------""" + cls.runModule( + "v.unpack", + input="data/tiny_polygons.pack", + output=cls.tiny_polygons, + ) + # This test vector contains a few small areas. One of the areas + # is so tiny that the sign of the signed area size might be wrong + # due to fp precision limits. In this case, the area is + # erroneously identified as an island and not as an area + cls.runModule( + "v.build", + map=cls.tiny_polygons, + ) + @classmethod def tearDownClass(cls): """Remove created vector map and temporary files, then delete the temp region.""" gs.run_command("g.remove", type="vector", flags="f", name="test_3x3_map") + gs.run_command("g.remove", type="vector", flags="f", name=cls.tiny_polygons) cls.del_temp_region() def test_vbuild_output(self): @@ -87,6 +104,26 @@ def test_vbuild_output(self): ) self.assertMultiLineEqual(filtered_output, self.vbuild_output) + def test_areas_islands(self): + """Compare the v.info output to the expected output.""" + self.assertModuleKeyValue( + "v.info", + map=self.tiny_polygons, + flags="tg", + sep="=", + precision=0.1, + reference={ + "nodes": 159, + "points": 0, + "lines": 0, + "boundaries": 168, + "centroids": 5, + "areas": 10, + "islands": 1, + "primitives": 173, + }, + ) + if __name__ == "__main__": test() From b55a971f87657af37dfab7e79f1623be839bfdf9 Mon Sep 17 00:00:00 2001 From: Markus Metz Date: Tue, 5 May 2026 14:41:50 +0200 Subject: [PATCH 3/7] improve code comment --- lib/gis/area_poly2.c | 12 +++++++----- lib/vector/diglib/poly.c | 12 +++++++----- 2 files changed, 14 insertions(+), 10 deletions(-) diff --git a/lib/gis/area_poly2.c b/lib/gis/area_poly2.c index 1cdecd92e7a..1d8cd2316e4 100644 --- a/lib/gis/area_poly2.c +++ b/lib/gis/area_poly2.c @@ -27,11 +27,13 @@ double G_planimetric_polygon_area(const double *x, const double *y, int n) double x1, y1, x2, y2, xshift, yshift; double area; - /* get a reference point - * do not overdo it, the first point would be good enough - * particularly for small polygons - * shift coordinates towards this reference point to make - * calculation of the signed area more robust */ + /* Get a reference point close to the polygon as new origin. + * Do not overdo it, the first point would be good enough, + * particularly for small polygons. + * Shift coordinates towards this reference point to make + * calculation of the signed area more robust by increasing the + * accuracy for given fp precision limits. + * See also dig_find_area_poly() in lib/vector/diglib/poly.c */ xshift = (x[0] + x[n / 2]) / 2.; yshift = (y[0] + y[n / 2]) / 2.; diff --git a/lib/vector/diglib/poly.c b/lib/vector/diglib/poly.c index c13ed8c84c5..fc0cde3eed1 100644 --- a/lib/vector/diglib/poly.c +++ b/lib/vector/diglib/poly.c @@ -103,11 +103,13 @@ int dig_find_area_poly(struct line_pnts *Points, double *totalarea) x = Points->x; y = Points->y; - /* get a reference point - * do not overdo it, the first point would be good enough - * particularly for small polygons - * shift coordinates towards this reference point to make - * calculation of the signed area more robust */ + /* Get a reference point close to the polygon as new origin. + * Do not overdo it, the first point would be good enough, + * particularly for small polygons. + * Shift coordinates towards this reference point to make + * calculation of the signed area more robust by increasing the + * accuracy for given fp precision limits. + * See also G_planimetric_polygon_area() in lib/gis/area_poly2.c */ mid = Points->n_points / 2; xshift = (x[0] + x[mid]) / 2.; yshift = (y[0] + y[mid]) / 2.; From b11723815357254734ad993886a52280611d0f5d Mon Sep 17 00:00:00 2001 From: Markus Metz Date: Tue, 5 May 2026 19:29:08 +0200 Subject: [PATCH 4/7] example to correctly compare fp numbers --- vector/v.to.db/testsuite/test_v_to_db.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/vector/v.to.db/testsuite/test_v_to_db.py b/vector/v.to.db/testsuite/test_v_to_db.py index 723af416de7..90add1863da 100644 --- a/vector/v.to.db/testsuite/test_v_to_db.py +++ b/vector/v.to.db/testsuite/test_v_to_db.py @@ -1,4 +1,5 @@ import json +import math from itertools import zip_longest from grass.gunittest.case import TestCase @@ -16,7 +17,12 @@ def _assert_dict_almost_equal(self, d1, d2): self.assertEqual(set(d1.keys()), set(d2.keys())) for k1 in d1: if isinstance(d1[k1], float): - self.assertAlmostEqual(d1[k1], d2[k1], places=6) + m_obs, e_obs = math.frexp(d1[k1]) + m_ref, e_ref = math.frexp(d2[k1]) + self.assertEqual(e_obs, e_ref) + # confident in the double fp calculations? + # use 12 instead of 7 decimal places + self.assertAlmostEqual(m_obs, m_ref, places=12) else: self.assertEqual(d1[k1], d2[k1]) From a0e168566b5f96e4018bdd1370e0b28f6209e491 Mon Sep 17 00:00:00 2001 From: Markus Metz Date: Wed, 6 May 2026 09:56:46 +0200 Subject: [PATCH 5/7] adjust code comments --- lib/gis/area_poly2.c | 8 +++++--- lib/vector/diglib/poly.c | 8 +++++--- vector/v.to.db/testsuite/test_v_to_db.py | 11 ++++++++--- 3 files changed, 18 insertions(+), 9 deletions(-) diff --git a/lib/gis/area_poly2.c b/lib/gis/area_poly2.c index 1d8cd2316e4..d4dd4b6c709 100644 --- a/lib/gis/area_poly2.c +++ b/lib/gis/area_poly2.c @@ -28,12 +28,14 @@ double G_planimetric_polygon_area(const double *x, const double *y, int n) double area; /* Get a reference point close to the polygon as new origin. - * Do not overdo it, the first point would be good enough, - * particularly for small polygons. + * The first point would be good enough, particularly for small + * polygons. The average of the first and the mid-point is a fast + * approximation of a reference point with reduced distance to the + * ring's vertices, also for larger rings. * Shift coordinates towards this reference point to make * calculation of the signed area more robust by increasing the * accuracy for given fp precision limits. - * See also dig_find_area_poly() in lib/vector/diglib/poly.c */ + * Keep in sync with dig_find_area_poly() in lib/vector/diglib/poly.c */ xshift = (x[0] + x[n / 2]) / 2.; yshift = (y[0] + y[n / 2]) / 2.; diff --git a/lib/vector/diglib/poly.c b/lib/vector/diglib/poly.c index fc0cde3eed1..3e4cd7a8259 100644 --- a/lib/vector/diglib/poly.c +++ b/lib/vector/diglib/poly.c @@ -104,12 +104,14 @@ int dig_find_area_poly(struct line_pnts *Points, double *totalarea) y = Points->y; /* Get a reference point close to the polygon as new origin. - * Do not overdo it, the first point would be good enough, - * particularly for small polygons. + * The first point would be good enough, particularly for small + * polygons. The average of the first and the mid-point is a fast + * approximation of a reference point with reduced distance to the + * ring's vertices, also for larger rings. * Shift coordinates towards this reference point to make * calculation of the signed area more robust by increasing the * accuracy for given fp precision limits. - * See also G_planimetric_polygon_area() in lib/gis/area_poly2.c */ + * Keep in sync with G_planimetric_polygon_area() in lib/gis/area_poly2.c */ mid = Points->n_points / 2; xshift = (x[0] + x[mid]) / 2.; yshift = (y[0] + y[mid]) / 2.; diff --git a/vector/v.to.db/testsuite/test_v_to_db.py b/vector/v.to.db/testsuite/test_v_to_db.py index 90add1863da..aa03070574b 100644 --- a/vector/v.to.db/testsuite/test_v_to_db.py +++ b/vector/v.to.db/testsuite/test_v_to_db.py @@ -20,9 +20,14 @@ def _assert_dict_almost_equal(self, d1, d2): m_obs, e_obs = math.frexp(d1[k1]) m_ref, e_ref = math.frexp(d2[k1]) self.assertEqual(e_obs, e_ref) - # confident in the double fp calculations? - # use 12 instead of 7 decimal places - self.assertAlmostEqual(m_obs, m_ref, places=12) + # use 12 instead of default 7 decimal places + # for reliable double precision fp results + self.assertAlmostEqual( + m_obs, + m_ref, + places=12, + msg=f"Comparing mantissas of {d1[k1]} and {d1[k1]}", + ) else: self.assertEqual(d1[k1], d2[k1]) From 5b2d49724c8d7cc03641b0dcca6b1cbf01e6ff54 Mon Sep 17 00:00:00 2001 From: Markus Metz Date: Wed, 6 May 2026 16:21:05 +0200 Subject: [PATCH 6/7] justify xshift --- lib/gis/area_poly2.c | 13 +++++++++++++ lib/vector/diglib/poly.c | 13 +++++++++++++ 2 files changed, 26 insertions(+) diff --git a/lib/gis/area_poly2.c b/lib/gis/area_poly2.c index d4dd4b6c709..5afd5cd82c1 100644 --- a/lib/gis/area_poly2.c +++ b/lib/gis/area_poly2.c @@ -32,10 +32,23 @@ double G_planimetric_polygon_area(const double *x, const double *y, int n) * polygons. The average of the first and the mid-point is a fast * approximation of a reference point with reduced distance to the * ring's vertices, also for larger rings. + * * Shift coordinates towards this reference point to make * calculation of the signed area more robust by increasing the * accuracy for given fp precision limits. + * + * Considering the basic formula + * area += (x2 - x1) * (y2 + 1)) + * the shift is in theory only needed for addition, not subtraction. + * Currently, x coordinates are subtracted and y coordinates are added. + * The reverse, adding x and subtracting y, is also correct, + * it simply reverts the sign. If for some reason the formula would + * be changed in the future by adding x and subtracting y, + * the current shift for both x and y would guarantee that + * the results stay correct (with reversed sign). + * * Keep in sync with dig_find_area_poly() in lib/vector/diglib/poly.c */ + xshift = (x[0] + x[n / 2]) / 2.; yshift = (y[0] + y[n / 2]) / 2.; diff --git a/lib/vector/diglib/poly.c b/lib/vector/diglib/poly.c index 3e4cd7a8259..729f3cd0a22 100644 --- a/lib/vector/diglib/poly.c +++ b/lib/vector/diglib/poly.c @@ -108,10 +108,23 @@ int dig_find_area_poly(struct line_pnts *Points, double *totalarea) * polygons. The average of the first and the mid-point is a fast * approximation of a reference point with reduced distance to the * ring's vertices, also for larger rings. + * * Shift coordinates towards this reference point to make * calculation of the signed area more robust by increasing the * accuracy for given fp precision limits. + * + * Considering the basic formula + * area += (x2 - x1) * (y2 + 1)) + * the shift is in theory only needed for addition, not subtraction. + * Currently, x coordinates are subtracted and y coordinates are added. + * The reverse, adding x and subtracting y, is also correct, + * it simply reverts the sign. If for some reason the formula would + * be changed in the future by adding x and subtracting y, + * the current shift for both x and y would guarantee that + * the results stay correct (with reversed sign). + * * Keep in sync with G_planimetric_polygon_area() in lib/gis/area_poly2.c */ + mid = Points->n_points / 2; xshift = (x[0] + x[mid]) / 2.; yshift = (y[0] + y[mid]) / 2.; From 30cea2970a33229cc6bf13e5ddf5c4d795d4f7d9 Mon Sep 17 00:00:00 2001 From: Markus Metz Date: Wed, 6 May 2026 16:59:05 +0200 Subject: [PATCH 7/7] fix comments, fix test message --- lib/gis/area_poly2.c | 11 +++-------- lib/vector/diglib/poly.c | 11 +++-------- vector/v.to.db/testsuite/test_v_to_db.py | 2 +- 3 files changed, 7 insertions(+), 17 deletions(-) diff --git a/lib/gis/area_poly2.c b/lib/gis/area_poly2.c index 5afd5cd82c1..cb9a24c341e 100644 --- a/lib/gis/area_poly2.c +++ b/lib/gis/area_poly2.c @@ -38,14 +38,9 @@ double G_planimetric_polygon_area(const double *x, const double *y, int n) * accuracy for given fp precision limits. * * Considering the basic formula - * area += (x2 - x1) * (y2 + 1)) - * the shift is in theory only needed for addition, not subtraction. - * Currently, x coordinates are subtracted and y coordinates are added. - * The reverse, adding x and subtracting y, is also correct, - * it simply reverts the sign. If for some reason the formula would - * be changed in the future by adding x and subtracting y, - * the current shift for both x and y would guarantee that - * the results stay correct (with reversed sign). + * area += (x2 - x1) * (y2 + y1) + * the shift is in theory only needed for addition, not subtraction, + * but does no harm and is kept for symmetry treating x and y coords. * * Keep in sync with dig_find_area_poly() in lib/vector/diglib/poly.c */ diff --git a/lib/vector/diglib/poly.c b/lib/vector/diglib/poly.c index 729f3cd0a22..ddd47a4e7e4 100644 --- a/lib/vector/diglib/poly.c +++ b/lib/vector/diglib/poly.c @@ -114,14 +114,9 @@ int dig_find_area_poly(struct line_pnts *Points, double *totalarea) * accuracy for given fp precision limits. * * Considering the basic formula - * area += (x2 - x1) * (y2 + 1)) - * the shift is in theory only needed for addition, not subtraction. - * Currently, x coordinates are subtracted and y coordinates are added. - * The reverse, adding x and subtracting y, is also correct, - * it simply reverts the sign. If for some reason the formula would - * be changed in the future by adding x and subtracting y, - * the current shift for both x and y would guarantee that - * the results stay correct (with reversed sign). + * area += (x2 - x1) * (y2 + y1) + * the shift is in theory only needed for addition, not subtraction, + * but does no harm and is kept for symmetry treating x and y coords. * * Keep in sync with G_planimetric_polygon_area() in lib/gis/area_poly2.c */ diff --git a/vector/v.to.db/testsuite/test_v_to_db.py b/vector/v.to.db/testsuite/test_v_to_db.py index aa03070574b..8516eaaa350 100644 --- a/vector/v.to.db/testsuite/test_v_to_db.py +++ b/vector/v.to.db/testsuite/test_v_to_db.py @@ -26,7 +26,7 @@ def _assert_dict_almost_equal(self, d1, d2): m_obs, m_ref, places=12, - msg=f"Comparing mantissas of {d1[k1]} and {d1[k1]}", + msg=f"Comparing mantissas of {d1[k1]} and {d2[k1]}", ) else: self.assertEqual(d1[k1], d2[k1])