From 366b782152c3897b2786b4d67a2c4a4b9faa3f89 Mon Sep 17 00:00:00 2001 From: rCarto Date: Mon, 27 May 2024 17:17:24 +0200 Subject: [PATCH] feat: add a smooth parameter to osrmIso*() functions, see #129 --- DESCRIPTION | 6 ++- R/osrmIsochrone.R | 79 ++++++++++++++++++++++++++------ R/osrmIsodistance.R | 55 +++++++++++++++++++++- R/utils.R | 13 ++++-- inst/tinytest/fill_grid_out.rds | Bin 7441 -> 9140 bytes man/osrmIsochrone.Rd | 9 ++++ man/osrmIsodistance.Rd | 9 ++++ 7 files changed, 150 insertions(+), 21 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index b932381..8b02589 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -25,7 +25,11 @@ Imports: sf Depends: R (>= 3.5.0) -Suggests: mapsf, tinytest, covr +Suggests: + mapsf, + tinytest, + covr, + terra URL: https://github.com/riatelab/osrm BugReports: https://github.com/riatelab/osrm/issues Encoding: UTF-8 diff --git a/R/osrmIsochrone.R b/R/osrmIsochrone.R index c3bb69c..a72ef61 100644 --- a/R/osrmIsochrone.R +++ b/R/osrmIsochrone.R @@ -19,6 +19,11 @@ #' grid, the total number of points will be res*res. Increase res to obtain more #' detailed isochrones. #' @param returnclass deprecated. +#' @param smooth if TRUE a moving window with a gaussian blur is applied to +#' durations. This option may be usefull to remove small patches of hard to +#' reach areas. The computed isochrones are less precise but better looking. +#' @param k size (sigma) of the gaussian moving window. A reasonable value is +#' used by default. #' @param osrm.server the base URL of the routing server. #' getOption("osrm.server") by default. #' @param osrm.profile the routing profile to use, e.g. "car", "bike" or "foot" @@ -64,27 +69,30 @@ #' } #' } osrmIsochrone <- function(loc, breaks = seq(from = 0, to = 60, length.out = 7), - exclude, res = 30, returnclass, + exclude, res = 30, smooth = FALSE, k, + returnclass, osrm.server = getOption("osrm.server"), osrm.profile = getOption("osrm.profile")) { opt <- options(error = NULL) on.exit(options(opt), add = TRUE) - + if (!missing(returnclass)) { warning('"returnclass" is deprecated.', call. = FALSE) } - + # input management loc <- input_route(x = loc, id = "loc", single = TRUE) oprj <- loc$oprj loc <- st_as_sf(data.frame(lon = loc$lon, lat = loc$lat), - coords = c("lon", "lat"), crs = 4326 + coords = c("lon", "lat"), crs = 4326 ) loc <- st_transform(loc, "epsg:3857") - + # max distance management to see how far to extend the grid to get measures breaks <- unique(sort(breaks)) tmax <- max(breaks) + + # gentle sleeptime & param for demo server if (osrm.profile %in% c("foot", "walk")) { speed <- 10 * 1000 / 60 } @@ -95,8 +103,8 @@ osrmIsochrone <- function(loc, breaks = seq(from = 0, to = 60, length.out = 7), speed <- 120 * 1000 / 60 } dmax <- tmax * speed - - + + # gentle sleeptime & param for demo server if (osrm.server != "https://routing.openstreetmap.de/") { sleeptime <- 0 @@ -105,7 +113,7 @@ osrmIsochrone <- function(loc, breaks = seq(from = 0, to = 60, length.out = 7), sleeptime <- 1 deco <- 75 } - + # create a grid to obtain measures sgrid <- rgrid(loc = loc, dmax = dmax, res = res) # slice the grid to make several API calls @@ -140,19 +148,20 @@ osrmIsochrone <- function(loc, breaks = seq(from = 0, to = 60, length.out = 7), listDur[[ltot]] <- dmat$durations listDest[[ltot]] <- dmat$destinations } - + measure <- do.call(c, listDur) destinations <- do.call(rbind, listDest) # for testing purpose # return(list(destinations = destinations, measure = measure, # sgrid = sgrid, res = res, tmax = tmax)) - + # assign values to the grid sgrid <- fill_grid( destinations = destinations, measure = measure, sgrid = sgrid, res = res, tmax = tmax ) - if (min(sgrid$measure) >= tmax + 1) { + + if (min(sgrid$measure, na.rm = TRUE) > tmax) { warning( paste0( "An empty object is returned. ", @@ -169,19 +178,63 @@ osrmIsochrone <- function(loc, breaks = seq(from = 0, to = 60, length.out = 7), ) return(empty_res) } + + if (isFALSE(smooth)) { + # All values not within breaks are set to tmax+1 + sgrid[is.na(sgrid$measure), "measure"] <- tmax + 1 + sgrid[is.nan(sgrid$measure), "measure"] <- tmax + 1 + sgrid[is.infinite(sgrid$measure), "measure"] <- tmax + 1 + sgrid[sgrid$measure > tmax, "measure"] <- tmax + 1 + } else { + if (!requireNamespace("terra", quietly = TRUE)) { + stop(paste0( + "'terra' package is needed for this function to work.", + "Please install it." + ), call. = FALSE) + } + r <- terra::rast(sgrid[, c("COORDX", "COORDY", "measure"), drop = TRUE], + crs = "epsg:3857") + if (missing(k)) { + k <- terra::res(r)[1] / 2 + } + mat <- terra::focalMat(x = r, d = k, type = "Gauss") + + # test for invalid focal matrix + if (sum(dim(mat)) < 6){ + warning( + paste0( + "An empty object is returned. ", + "Select a larger value for 'k'." + ), + call. = FALSE + ) + empty_res <- st_sf( + crs = ifelse(is.na(oprj), 4326, oprj), + id = integer(), + isomin = numeric(), + isomax = numeric(), + geometry = st_sfc() + ) + return(empty_res) + } + sgrid <- terra::focal(x = r, w = mat, fun = mean, na.rm = TRUE) + sgrid[is.na(sgrid)] <- tmax + 1 + } + + # computes isopolygones iso <- mapiso(x = sgrid, breaks = breaks, var = "measure") # get rid of out of breaks polys iso <- iso[-nrow(iso), ] # fisrt line always start at 0 iso[1, "isomin"] <- 0 - + # proj mgmnt if (!is.na(oprj)) { iso <- st_transform(x = iso, oprj) } else { iso <- st_transform(x = iso, 4326) } - + return(iso) } diff --git a/R/osrmIsodistance.R b/R/osrmIsodistance.R index bb3b897..30e6712 100644 --- a/R/osrmIsodistance.R +++ b/R/osrmIsodistance.R @@ -19,6 +19,11 @@ #' square grid, the total number of points will be res*res. Increase res to #' obtain more detailed isodistances. #' @param returnclass deprecated. +#' @param smooth if TRUE a moving window with a gaussian blur is applied to +#' distances. This option may be usefull to remove small patches of hard to +#' reach areas. The computed isodistances are less precise but better looking. +#' @param k size (sigma) of the gaussian moving window. A reasonable value is +#' used by default. #' @param osrm.server the base URL of the routing server. #' getOption("osrm.server") by default. #' @param osrm.profile the routing profile to use, e.g. "car", "bike" or "foot" @@ -62,7 +67,8 @@ #' } #' } osrmIsodistance <- function(loc, breaks = seq(from = 0, to = 10000, length.out = 4), - exclude, res = 30, returnclass, + exclude, res = 30, smooth = FALSE, k, + returnclass, osrm.server = getOption("osrm.server"), osrm.profile = getOption("osrm.profile")) { opt <- options(error = NULL) @@ -138,7 +144,9 @@ osrmIsodistance <- function(loc, breaks = seq(from = 0, to = 10000, length.out = destinations = destinations, measure = measure, sgrid = sgrid, res = res, tmax = tmax ) - if (min(sgrid$measure) >= tmax + 1) { + + + if (min(sgrid$measure, na.rm = TRUE) >= tmax) { warning( paste0( "An empty object is returned. ", @@ -155,6 +163,49 @@ osrmIsodistance <- function(loc, breaks = seq(from = 0, to = 10000, length.out = ) return(empty_res) } + + if (isFALSE(smooth)) { + # All values not within breaks are set to tmax+1 + sgrid[is.na(sgrid$measure), "measure"] <- tmax + 1 + sgrid[is.nan(sgrid$measure), "measure"] <- tmax + 1 + sgrid[is.infinite(sgrid$measure), "measure"] <- tmax + 1 + sgrid[sgrid$measure > tmax, "measure"] <- tmax + 1 + } else { + if (!requireNamespace("terra", quietly = TRUE)) { + stop(paste0( + "'terra' package is needed for this function to work.", + "Please install it." + ), call. = FALSE) + } + r <- terra::rast(sgrid[, c("COORDX", "COORDY", "measure"), drop = TRUE], + crs = "epsg:3857") + if (missing(k)) { + k <- terra::res(r)[1] / 2 + } + mat <- terra::focalMat(x = r, d = k, type = "Gauss") + + # test for invalid focal matrix + if (sum(dim(mat)) < 6){ + warning( + paste0( + "An empty object is returned. ", + "Select a larger value for 'k'." + ), + call. = FALSE + ) + empty_res <- st_sf( + crs = ifelse(is.na(oprj), 4326, oprj), + id = integer(), + isomin = numeric(), + isomax = numeric(), + geometry = st_sfc() + ) + return(empty_res) + } + sgrid <- terra::focal(x = r, w = mat, fun = mean, na.rm = TRUE) + sgrid[is.na(sgrid)] <- tmax + 1 + } + # computes isopolygones iso <- mapiso(x = sgrid, breaks = breaks, var = "measure") # get rid of out of breaks polys diff --git a/R/utils.R b/R/utils.R index 16e50f3..841fca5 100644 --- a/R/utils.R +++ b/R/utils.R @@ -372,16 +372,19 @@ fill_grid <- function(destinations, measure, sgrid, res, tmax) { n = c(res, res) ) ag_pt <- function(x) { - if (length(x > 0)) { + if (length(x) > 0) { min(rpt[["measure"]][x], na.rm = TRUE) } else { - tmax + 1 + NA } } inter <- st_intersects(xx, rpt) sgrid$measure <- unlist(lapply(inter, ag_pt)) - sgrid[is.infinite(sgrid$measure), "measure"] <- tmax + 1 - sgrid[is.nan(sgrid$measure), "measure"] <- tmax + 1 - sgrid[sgrid$measure > tmax, "measure"] <- tmax + 1 + sgrid[is.infinite(sgrid$measure), "measure"] <- NA + sgrid[is.nan(sgrid$measure), "measure"] <- NA sgrid } + + + + diff --git a/inst/tinytest/fill_grid_out.rds b/inst/tinytest/fill_grid_out.rds index bbae8e02dffcc68bf67a3839082975c5f8ea7987..5b62bfcfbe55ab87169174912b15fbe30081f966 100644 GIT binary patch literal 9140 zcmchcc~p~E`p2E0?bw#q*t#Ghq^%28ic}F1S>9T;SP`f)C9);cT4gr{gplv1i{}@UPE8RS7*YT26=JB7@zCsdhegx4q*Ig(ny9&Rw_f-Zcr?t?O=Lc&_2Q zz_4XSxu2Y-e6&t`aP7K%&*4hx3Wu*<0yA$!9(a0a^FzCcPC?oIHImUTCId*gG3-R}2LQ+>^ z{xIW`DqvXBX$<3F55t=|4K=*;s`g<MM2B=&jA$W^%vGtQ{Wo=Q55 z!5r36IF3WA;hj~DC9027F=A?AEFH3DZ!z@=@@+i5`Z|b3;FnMU%CMOF z1uPq)$NqvYKo-3OK1XZyK;XZTss{qle_^Ku(+6I#{iMfULzf_#FM-{tP7j3q8{b+m z{oobb13j2ybMHg?zv8)x-`Lm^!ZcmCck zn>*v;P%==#={ZL{O7oBwIY8YrN87wS_yw8(2k_Ew(Y_Kqr}-@$z+z4Nf5Pn?@mrdm zG;$@>GP9zs;WPdzjmt{#$NvWhFjwR9pD0@ip3(dURp*FDXs*&o2dHxZ72k-?`N+KDVykiHiDWO;>DGN8V}a54B)Kwlxj3erKrMawIsl<-bg# z2h`UX_6BY_vYEHN1GfSSLXIg&p@JgID9u{{IGjKz*Ho zHNfRi0OGGO9S~$&ULoPi!Bdvsg9z)d-<G(rABCts8>X) zfRME45%9U~Pg=|c@Nd?I6W9;Y1?IU+9l&1OPAw)EykPxl0=pLd(7b4_QsyIIx9uY> zCI|exb=L&;V|0mm?otP^-}bE*a}m5^eK3JtkA7y(oU0V@2>8PGzBc;;c+$FN0{apA zmU-?{2e8j}K%1QlUb0@Fz^+3-F$c_5YJUWLY1^#L&H>L^4^Lp%pexODmpXu8Te3F$ zB3NvdPGDWo?dJBmN-mFpowm2NAs4_OtxG3ZAEIxV=Pq>sU)jFWhU9{|){7IYwdh*2 z%Uq?hN5CH225m?Vc*Yt!!TK0oW}dS|_!teG08d-pbI{Uw+YW!1PxISCnzc0f-IU2L zFa-$_t#-?|+3xmT7G#u0Sro)`)>_b(-lWHX#-D`xv9L{uRDPM3mRb=-F`Vy&Zux|6 zLi$5^wcNpC9a*mHm6r&8 zSR{WnO=52Tox)gitNf(UE78Bs+?~_ZUrjkOybxo}o}Dq+)9#XKrccu3X*Gy;mZX+N z!(uj{KjoTQmSiQ+2{uicJebRXnEmdkWHippD#=wv4hc)r`+qcHX(qz#OqbK3s27$PpM;7*ok39OqZ$Gjzvx~*9F1#mqmYDj))gO>)QPx> zt|Q5UxkwdmMOavGNqT#3zbKy;GhG*gAGVHReLYDdWsD59*O!c~v$K4H+O=f%i`6w= zi}C32XdDThc;U&`&;C$BT7a7VoY~YrjdT{W_^4ra2G^b<@{GTZ#$gia^1)pKfHV>V zDh)Z@hTalm6M&6lNU5o1esz@=@4nI4N)UbN^q?Od4aMO|RE_rrR@va7TO-|K5Gz&; zxD_8$`=)c6Qf2JULNj}*`dOnRj;mz0Llk zFShc7-@m3QJ-6gMJN;xVZRgrqll8MEVQ)W~=?t7TiP^GMz*u@svyM4y@_(mdX;(j) zHTi1Ro3R$#Lk>8Tb(xS z>i;8#FI-o@dZk*)@f!+o=3>f?cY^{tA0Ph%hBJ@fdCf|<>}4*leq7h` zn)6Wx1>CZ|xn{#__A(b&->UlI)#uYq>Cd?ZZU~>hkP*HwoWGFev*s~mn9ze+m#$9k zbtwEJhW8`JU%gU+;pUoeP?v=H%MbMM(52T@{9GE!$*}FOZbnCWlY0lz0_waaV-;xq zT>j+ryg4Uc!Q;yLAvfl$jZm7UJ*7lj@Vp(UipaZr=kAslug=`L6vncv%U@GX;lv$% zS?kK!?&kUP&xM_6o5nErdA~Y==W@mV*V|M-fKelDkFw{jvH0<<6vE{0c}tk+!^G9i zw%U2CKrX1&7UC?SugR&obnpM#sQqfEU-Cjm=j~rgoEkS|)g$bAyY@y*F?J9|H9l_< zT*Tq)6jz~;c`LywBF{g?G=Dds5Hau7IYNQ+<>2vH(bb-qta-VVY)9(~4x}|`wC?r1 z$0^gaQj#Sbl>d|u%oNmjwdkDO+9P}-#p0V@(@iz~K|Ye?+t*Vqkf~??k6$f6cUf2B zB4Wo!)pA~z7fm+mq*UWHBSk5pXW9b|bOK$fummVf$2rHHhy>&D!1fy1K5F~(Qj#bH z^^K4;B*C09F;M63=K*BrnXnlYi@}&x+pcVAsKXK~>YOVfkvKfg8z)E^yG@T_a8#pB zhlTtxA7}3Uja3oi#+1R4()X0B^1W_{iyz;K4~iR{I^Y5LG^j@S=;?*r^eLcuCMvM> z`VmcEUO?9BrWz@e=fjk2^UGH`xw$2cRka7Z;qQI7?PhZV+?(`Gi{|pher}MuhSu6y zUn#mBe@Jxa=!<)4mAER)XDLSu#l|TyKVwqy&3o&poN0++io=7`rzB|^8m3Vf9aH3^ zU?r)h;{C-V!%+qP(GA;kJ%A9R7=1Y1&d@YkvR8Uhx4GiB-exZ6H8=PMNjRoopO^IC z*eB9lwiuThl6vy0vt+dfl~a{Vq}W!NDGf`FO^mA1&JIK$Gn(|ldkt3}N_%2_A{CK^ zKCV?sF2s#+PUWQiUp62u8x&bTbV451k}7PHF&)_< zsQG$X+v?)TC&Y)RN?fvdU+doW2^#}t-pAX2=$DlcLs08vz3C`(k3gSY*eRfXCo7BV zIvvl|BzOSmFL649!gx~=S{Wz3rcgRQB`$=MA`JIK`YCvV(B`;Oz!@P7Buuf5?5oj{ zB_W|}RI3^*$8X!x^8|&qgZ07p#c#x{G0Z-v$qWGnrlH=A^5n)=-yXNJ2*ESr?E*pj zku|!AcAAJEv>5NMp{TaO9j=Y*Wc_am9wxhW@LHR_SwThP#H;v=Y%Wlo#&ofW3_^t< z4)1s+(|c7yYQLI1`cPx^rAoqT(2WXqMYuZ??VBPmQy!!Dkf9c(OuXM;bpJ~{@jK`9 zy|wq~p}k*b$l_1pCQ_rd7DXX>+XxY#;aVWvx7|6@)3ZI4(PJ$fqR`5@9m+mswua9? zwL=+St@~0oPTrunPvTX5D_NKNk7%hORv=oZ+Sq1KPUvkA=~j^?*eJ||Ym@Uh*`Ej9 zjwx0O4z-J7z!BIy2)Ocvn*9!>n&91g%2eJWSr7}}m-!gjkD7$OKOh17D30;i9 zZuvGuW`0-PnmleO?Y`b9@e|Z7?lkVuf7u5uRj| zMZR0cL>TAEpMHxHqz2(rO7TK#;-@yzc^>0Z8Q*)Qj_QE$_7?TVo^R}>B(6}BI(gRbIt+90c zfJ|>5r)|f!BIBV0o+92L!m%!gA?jJ3(tuttZBU0uaK3;KV7WK2v>8MEm2yhjn3o6_ zn;3U~6!F|GzK0?_1j4~$7NIl^k?e_VS7KvmIk?~+`u!c0Iw6UEU#{8_!K3&gr-150 zI$=l|)h=O20Wzca%T#%-XOdr1GQ+_{8s>|XO|oVRu#xHIOcd3K)+CLGPL;W^%OF({ zUS3S{6Q@0=^M2i>MyrW45-{IlEhjSvLP@hv=?YqI66_>TZI|B>dW&j=dmz@#=myE6 z&4lcL8m5C9-Y8_7+|&zOkxC=~1+ADRa)6a~{y+~A&+>{=x6*syUk|F!GNC=SxNqPeEZH|oErIykif}OhUZa3KM)i#X^{`cIRbG`qxW;3{M#~S#GJ&X~ zBiuJe;78>uxRJOEtE9uavWKq1ZR{0pVUMqyp%a<}V za0R=EPgrJ()-_2UT}SOjYpO(5I!Qzos8feVV>knv#+NQi-m%1l?1j=Z${xcZ9Yxd3 zj!S@`E47W-;h}DMn81Z$y^3k{;_(ZlLvLh5W%eSTAOL!H&otAg_SK(wZdzN{3J0?n zWoNQs+MO3O%S=m!SlC`bx&cA%MC5vWWVI}D;{(@zB7eOCZu+|j)WGqLt#nLQlZc8o z`ooMv#U!@9u@XU&G}NuxWn}MOJjpRZCO0$$}z3-C9{t9 za!qfm7p6oJCAc&46Nobg9^T{&N%xNhfMF_=Z{k&!Fz{GJ1LN3)c^fI2=DQf=?ohjR zL#p*dT`PJ+l1EGOMSiPU3ys{xv@xMPo_XDC+O=MuJ^E?)&Go`Sd4N3WmNTnR=%&^? zZt-!JhG50_s3uAK;@%Oh*n?D(7{@X%Xc7@Sixc9erxdXjNc#Q&swoBdrNsPBk7}?< zvN?)-Sbo|*tDdFWAqw+?4?lK170(aur`p85-C|FqzrLHF?f literal 7441 zcmchZ2~d;Q7RTFCp|m0{t(48ESa=o@1SCY3v|@`^k>?7OKv?w#+Rd(Qvd zdqJ5Wz4xz2hgXRHaWDZ*DYWW8_5Igpu&lD5q8)B{9msK}Lbh8_Uq%orxy>=p%k&4_mHuL(r%QR{fX-3>j4d4( zm^9%GIw?H$CsO@y&9%qf)iM-;p%@n0(mom5&GYp+cC7ZdFiFQ}yYKYhaIf-H#r}TcqC9n??Dx|GPVkMz~f!JwE?^c4nAY{1h?&J8W_f~9S<*a5>@ zP#Cleb4wK{MGXi}lXKsKyHvpxgp#vWZyR6A;3k|8wi$S2S8!KE!7)XC-roHbiv`lWqoAnP5p#Z7xaONV~1kCC3S5> zPOuM6l4@K%sh}@_8_FZ1cYwu&XM6S)&=%-k%mw8~vUQ?A8A>KRsoB>+UqHXZ+)(~a zc0}|iL+AvhJ^L!i0{SiH7v<-&U7|l3DD5B%;1}gC(Yt^OvH}W}(IU$7VAW384E(Ir z7rhIpgl)hzUp|gHI=tZ1`a! zY&kuh7YBx z=@7fQ0q-)+tZm?GUGCi)NDmq% zgIg14EBk{+DcQ)ynZGa1&Ydwb;kSPJqwi2WcK|@SG4D(7%cQw>ZtiA{)l<^LYmOAN4v=-=yb{R6$rL?ng!|eiqD*9Djc_7IMMwUKn(Xp0rNajnB^@|J8V{|W#D0lz zpD&bnmy!J;;`V;)+LYlAoM&zBo@(Xon6S(pT4`*e*Sq8&3b#Emyy4W;Eo)I{chwFfA00e`pKT>89CBt!KvG8Yd z#xu6=A;;nKLB+z2x*INo_udL87CzZG?qzIgg5=pfPfWeIH003Ikh>wp!W^TeA+*Hs zWX0u-XZeXs){jpW3)8o5xZL;!E}Ymz-8XI%T(0nXQm0HWNcir$MN)XoHbbwhi6ZJ% zS4YS-tfi-@XJY@RJTt1tf8F_@r(?Mv@xhQ@QGQ~;a{YhJi@nw^%Xr*0GO@o?^D4V@ z<*WKbJp95>FBhJLe6uq$_{R7ZhVT!Q{WCM4RLw4v3l8SGg^gyF#m?m%`b)_?RdmyQ znZyqwV@fBMT!3r#y|!0vbLMik^1t6hy^{j^;5yArgS!lUD0r+qCVt=3_A3mn-xo~qS3^|t1M z;p$>-Zm5;Aw9JT5wpQhXb;6S$2xqh#T8`9RBl6JdoLEw=x^g~2tFiTL2*er8Yjs1< z3bT!a$x^hc+DVi*zxm)P(-8A3)tv1iG-bzrpM=pt{P1(M7Wi>54G1k^H zTOC-vC%Zo;DXl#aZK)PDmAw_@cuA)Z^lV}~>+1f8El?KJ^?36`6)SaE^ytLU)o?B# z@J(LtT4`_Lx}Mj9r0@d_CR6juQgz{&DTa9G z#NyBRTs{Iv5CowpWc6r+vMo^WUMGd&5w=CZydU)zbVpvu

b{=&^o zCzGG{N3u=h`K%-43JgQ)hKlDiqupY-XPoLAz>LV@rV{tt*PrE-Smm==hV^D{0gOGa z_D*Xr_2Fl;v(j|sB}Kkmc3dW$T9$3$MKNqJTO2^uBWn}JYX?YG7>>cq)X2x{n<{zd zhQ`)!ck16zZO*SU|M>T3TZ>XyWBUQPrpH#+8GVOyAGOfHYfN6#=7H8Ay^GZwy>*zcCdPmK@z*cjj+UO@jxSJ5`MFRVbU7!E$G~u{CU6d56)eEvQ55%Ao4t?=eIs}wvA@{@xqsWfB0nhTX8Pc znAv;0R(Xx4#C9k8e1zs(8nlBk!sEyeD~9g4!$#U90rRbyGlnk~EXd?}p`Ji4)+BLS zT?Hs^2}kpFHT0<@!1b&SL*7gh2;#Hh=P(xq70Fd_evQCt-VewW2wGG`;CZEhj>R}4G7&DTYV!YtFV8%cBMJ>Nvz;=T%NSJ;(~l@nI<5lIxr+qL$?V$ zf`%i1b#XFl(w`@#%q1p%WxII>y|*#W(D0kkqAjiIM5jK6k^VwjNU+;&vX$$nGud84 uuO>$bf?OfXmTd#)1Q_!sYvbICe8VG