From e0bfbccf0eba3e2fb0404a4a7c7ff421c1919888 Mon Sep 17 00:00:00 2001 From: Timon Date: Tue, 24 Feb 2026 13:07:55 +0100 Subject: [PATCH] Refactor area calculation functions and enhance error handling; add area conversion utility --- r_app/80_utils_agronomic_support.R | 23 +++++++--- r_app/80_utils_cane_supply.R | 43 ++++++++---------- r_app/80_utils_common.R | 43 +++++++----------- ..._CI_report_with_kpis_agronomic_support.Rmd | 16 +++---- r_app/parameters_project.R | 31 ++++++++++++- r_app/translations/translations.xlsx | Bin 44729 -> 44690 bytes 6 files changed, 89 insertions(+), 67 deletions(-) diff --git a/r_app/80_utils_agronomic_support.R b/r_app/80_utils_agronomic_support.R index 04f0e9a..9ee5cb2 100644 --- a/r_app/80_utils_agronomic_support.R +++ b/r_app/80_utils_agronomic_support.R @@ -356,8 +356,8 @@ calculate_tch_forecasted_kpi <- function(field_statistics, harvesting_data = NUL #' Uses FOUR_WEEK_TREND_* thresholds defined in 80_utils_common.R: #' - FOUR_WEEK_TREND_STRONG_GROWTH_MIN (0.3) #' - FOUR_WEEK_TREND_GROWTH_MIN (0.1) -#' - FOUR_WEEK_TREND_STABLE_THRESHOLD (0.1 and -0.1) -#' - FOUR_WEEK_TREND_WEAK_DECLINE_THRESHOLD (-0.3) +#' - FOUR_WEEK_TREND_STABLE_THRESHOLD (0.1; used as lower bound via -FOUR_WEEK_TREND_STABLE_THRESHOLD) +#' - FOUR_WEEK_TREND_WEAK_DECLINE_THRESHOLD (-0.1) #' - FOUR_WEEK_TREND_STRONG_DECLINE_MAX (-0.3) #' calculate_growth_decline_kpi <- function(ci_values_list) { @@ -632,14 +632,25 @@ create_field_detail_table <- function(field_boundaries_sf, all_kpis, current_wee # ============================================ # GROUP 0b: FIELD AREA (from geometry) # ============================================ + # Precompute unit label before tryCatch to avoid re-raising errors in error handler + unit_label <- get_area_unit_label(AREA_UNIT_PREFERENCE) + area_col_name <- paste0("Area_", unit_label) + # Calculate field area using unified function and preferred unit tryCatch({ field_areas <- calculate_area_from_geometry(field_boundaries_sf, unit = AREA_UNIT_PREFERENCE) - unit_label <- get_area_unit_label(AREA_UNIT_PREFERENCE) - area_col_name <- paste0("Area_", unit_label) + + # Validate that calculated areas match the number of fields + n_fields <- nrow(field_boundaries_sf) + if (length(field_areas) != n_fields) { + stop( + "Mismatch: calculate_area_from_geometry returned ", length(field_areas), + " areas but field_boundaries_sf has ", n_fields, " fields" + ) + } area_df <- data.frame( - field_idx = seq_along(field_areas), + field_idx = seq_len(n_fields), area_value = field_areas, stringsAsFactors = FALSE ) @@ -649,7 +660,7 @@ create_field_detail_table <- function(field_boundaries_sf, all_kpis, current_wee left_join(area_df, by = "field_idx") }, error = function(e) { message(paste("Warning: Could not calculate field areas:", e$message)) - result[[paste0("Area_", get_area_unit_label(AREA_UNIT_PREFERENCE))]] <<- NA_real_ + result[[area_col_name]] <<- NA_real_ }) # ============================================ diff --git a/r_app/80_utils_cane_supply.R b/r_app/80_utils_cane_supply.R index eaa904b..ca896a1 100644 --- a/r_app/80_utils_cane_supply.R +++ b/r_app/80_utils_cane_supply.R @@ -53,12 +53,22 @@ CI_CHANGE_INCREASE_THRESHOLD <- CI_CHANGE_RAPID_GROWTH_THRESHOLD # Weekly CI calculate_field_acreages <- function(field_boundaries_sf, unit = AREA_UNIT_PREFERENCE) { tryCatch({ # Get field identifier (handles pivot.geojson structure) - field_names <- field_boundaries_sf %>% - sf::st_drop_geometry() %>% - pull(any_of(c("field", "field_id", "Field_id", "name", "Name"))) - - if (length(field_names) == 0 || all(is.na(field_names))) { - field_names <- seq_len(nrow(field_boundaries_sf)) + # Use intersect() to find the first matching column, avoiding pull() ambiguity + field_names <- { + available_cols <- intersect( + colnames(field_boundaries_sf), + c("field", "field_id", "Field_id", "name", "Name") + ) + + if (length(available_cols) > 0) { + # Extract the first matching column + field_boundaries_sf %>% + sf::st_drop_geometry() %>% + pull(available_cols[1]) + } else { + # Fall back to sequential indices if no field name column exists + seq_len(nrow(field_boundaries_sf)) + } } # Use unified area calculation function @@ -70,15 +80,15 @@ calculate_field_acreages <- function(field_boundaries_sf, unit = AREA_UNIT_PREFE result_df <- data.frame( field = field_names, - area = areas, - stringsAsFactors = FALSE + warning(paste("Could not calculate areas from geometries -", e$message)) stringsAsFactors = FALSE ) colnames(result_df) <- c("field", col_name) # Aggregate by field to handle multi-row fields (e.g., sub_fields) + # Use bare lambda to preserve original column name (not list() which creates suffixes) result_df %>% group_by(field) %>% - summarise(across(all_of(col_name), list(~ sum(., na.rm = TRUE))), .groups = "drop") + summarise(across(all_of(col_name), ~ sum(., na.rm = TRUE)), .groups = "drop") }, error = function(e) { message(paste("Warning: Could not calculate areas from geometries -", e$message)) unit_label <- get_area_unit_label(unit) @@ -600,21 +610,6 @@ calculate_field_analysis_cane_supply <- function(setup, reports_dir ) - # cat("\n--- Per-field Results (first 10) ---\n") - # available_cols <- c("Field_id", "Acreage", "Age_week", "Mean_CI", "Four_week_trend", "Status_Alert", "Cloud_category") - # available_cols <- available_cols[available_cols %in% names(field_analysis_df)] - # if (length(available_cols) > 0) { - # print(head(field_analysis_df[, available_cols], 10)) - # } - - # ========== PHASE 10: CALCULATE FARM-LEVEL KPIS ========== - # farm_kpi_results <- calculate_farm_level_kpis( - # field_analysis_df, - # current_week, - # current_year, - # end_date - # ) - # For now, farm-level KPIs are not implemented in CANE_SUPPLY workflow farm_kpi_results <- NULL diff --git a/r_app/80_utils_common.R b/r_app/80_utils_common.R index e47bab0..f74b7d9 100644 --- a/r_app/80_utils_common.R +++ b/r_app/80_utils_common.R @@ -86,12 +86,12 @@ PHASE_DEFINITIONS <- data.frame( #' Calculate field area from geometry in specified unit #' -#' Unified function for calculating polygon area from sf or SpatVect geometries. +#' Unified function for calculating polygon area from sf or SpatVector geometries. #' Uses equal-area projection (EPSG:6933) for accurate calculations across all zones. #' -#' @param geometry sf or SpatVect object containing field polygons +#' @param geometry sf or SpatVector object containing field polygons #' If sf: must have geometry column (auto-detected) -#' If SpatVect: terra object with geometry +#' If SpatVector: terra object with geometry #' @param unit Character. Output unit: "hectare" (default) or "acre" #' #' @return Numeric vector of areas in specified unit @@ -118,7 +118,7 @@ PHASE_DEFINITIONS <- data.frame( #' areas_ha <- calculate_area_from_geometry(fields_sf, unit = "hectare") #' areas_ac <- calculate_area_from_geometry(fields_sf, unit = "acre") #' -#' # With SpatVect +#' # With SpatVector #' library(terra) #' fields_vect <- vect("pivot.geojson") #' areas_ha <- calculate_area_from_geometry(fields_vect, unit = "hectare") @@ -137,12 +137,12 @@ calculate_area_from_geometry <- function(geometry, unit = "hectare") { # Handle sf object geometry_proj <- sf::st_transform(geometry, 6933) areas_m2 <- as.numeric(sf::st_area(geometry_proj)) - } else if (inherits(geometry, "SpatVect")) { - # Handle terra SpatVect object + } else if (inherits(geometry, "SpatVector")) { + # Handle terra SpatVector object geometry_proj <- terra::project(geometry, "EPSG:6933") areas_m2 <- as.numeric(terra::expanse(geometry_proj)) } else { - stop("geometry must be an sf or terra SpatVect object. Got: ", + stop("geometry must be an sf or terra SpatVector object. Got: ", paste(class(geometry), collapse = ", ")) } @@ -1068,12 +1068,14 @@ calculate_field_statistics <- function(field_boundaries_sf, week_num, year, } # Guard: detect cloud-masked dates (CI == 0 indicates no-data) - # When any extracted value is 0, treat the entire date as cloud-masked - has_zeros <- any(extracted$CI == 0, na.rm = TRUE) + # Calculate proportion of zero values; only treat as cloud-masked if proportion exceeds threshold + prop_zero <- mean(extracted$CI == 0, na.rm = TRUE) + cloud_threshold <- 0.25 # 25% tolerance for zero pixels - if (has_zeros) { + if (prop_zero > cloud_threshold) { # Cloud-masked date: skip temporal analysis, set stats to NA - message(paste(" [CLOUD] Field", field_name, "- entire date is cloud-masked (CI==0)")) + message(paste(" [CLOUD] Field", field_name, "- date is cloud-masked (", + round(prop_zero * 100, 1), "% CI==0)")) results_list[[length(results_list) + 1]] <- data.frame( Field_id = field_name, @@ -1089,7 +1091,7 @@ calculate_field_statistics <- function(field_boundaries_sf, week_num, year, next } - ci_vals <- extracted$CI[!is.na(extracted$CI)] + ci_vals <- extracted$CI[!is.na(extracted$CI) & extracted$CI > 0] if (length(ci_vals) == 0) { next @@ -1847,21 +1849,8 @@ calculate_yield_prediction_kpi <- function(field_boundaries, harvesting_data, cu pred_rf_current_season$predicted_Tcha < yield_quartiles[3], na.rm = TRUE) lowest_25_count <- sum(pred_rf_current_season$predicted_Tcha < yield_quartiles[1], na.rm = TRUE) - # Calculate total area - if (!inherits(field_boundaries, "SpatVector")) { - field_boundaries_vect <- terra::vect(field_boundaries) - } else { - field_boundaries_vect <- field_boundaries - } - - # Handle both sf and SpatVector inputs for area calculation - if (inherits(field_boundaries, "sf")) { - field_boundaries_projected <- sf::st_transform(field_boundaries, "EPSG:6933") - field_areas <- sf::st_area(field_boundaries_projected) / 10000 # m² to hectares - } else { - field_boundaries_projected <- terra::project(field_boundaries_vect, "EPSG:6933") - field_areas <- terra::expanse(field_boundaries_projected) / 10000 - } + # Calculate total area using centralized function + field_areas <- calculate_area_from_geometry(field_boundaries, unit = "hectare") total_area <- sum(as.numeric(field_areas)) safe_log(paste("Total area:", round(total_area, 1), "hectares")) diff --git a/r_app/90_CI_report_with_kpis_agronomic_support.Rmd b/r_app/90_CI_report_with_kpis_agronomic_support.Rmd index b9ff3ab..0fc4b53 100644 --- a/r_app/90_CI_report_with_kpis_agronomic_support.Rmd +++ b/r_app/90_CI_report_with_kpis_agronomic_support.Rmd @@ -531,19 +531,20 @@ map_trend_to_arrow <- function(text_vec, include_text = FALSE) { } # Determine category and build output with translated labels - if (grepl("strong growth", text)) { + # Using word-boundary anchored patterns (perl=TRUE) to avoid substring mis-matches + if (grepl("\\bstrong growth\\b", text, perl = TRUE)) { arrow <- "↑↑" trans_key <- "Strong growth" - } else if (grepl("slight growth|weak growth|growth|increasing", text)) { + } else if (grepl("\\b(?:slight|weak) growth\\b|(?% merge_v(j = tr_key("KPI")) %>% autofit() @@ -1864,7 +1864,7 @@ metadata_info <- data.frame( format(Sys.time(), "%Y-%m-%d %H:%M:%S"), paste(tr_key("project"), toupper(project_dir)), paste(tr_key("week"), current_week, tr_key("of"), year), - ifelse(exists("AllPivots0"), nrow(AllPivots0 %>% filter(!is.na(field)) %>% group_by(field) %>% summarise()), tr_key("unknown")), + if (exists("AllPivots0")) { nrow(AllPivots0 %>% filter(!is.na(field)) %>% group_by(field) %>% summarise()) } else { tr_key("unknown") }, tr_key("next_wed") ) ) diff --git a/r_app/parameters_project.R b/r_app/parameters_project.R index d3af7bc..45d65a3 100644 --- a/r_app/parameters_project.R +++ b/r_app/parameters_project.R @@ -64,8 +64,35 @@ get_area_unit_label <- function(unit = AREA_UNIT_PREFERENCE) { } else if (unit_lower == "acre") { return("ac") } else { - warning("Unknown unit: ", unit, ". Defaulting to 'ha'") - return("ha") + warning("Unknown unit '", unit, "'. Falling back to AREA_UNIT_PREFERENCE ('", AREA_UNIT_PREFERENCE, "')") + return(get_area_unit_label(AREA_UNIT_PREFERENCE)) + } +} + +#' Get area conversion factor from hectares to requested unit +#' +#' Returns the numeric multiplier to convert 1 hectare into the requested unit. +#' Case-insensitive unit matching. On unknown unit, warns and defaults to 1 (hectare). +#' +#' @param unit Character. Unit preference ("hectare" or "acre") +#' @return Numeric. Conversion factor (1 for hectare, 2.47105 for acre) +#' +#' @details +#' Conversion factors: +#' - "hectare" → 1 (identity) +#' - "acre" → 2.47105 (1 hectare ≈ 2.47105 acres) +#' +#' +#' @export +get_area_conversion_factor <- function(unit = AREA_UNIT_PREFERENCE) { + unit_lower <- tolower(unit) + if (unit_lower == "hectare") { + return(1) + } else if (unit_lower == "acre") { + return(2.47105) # 1 hectare = 2.47105 acres (inverse of 0.404686) + } else { + warning("Unknown unit '", unit, "'. Falling back to AREA_UNIT_PREFERENCE ('", AREA_UNIT_PREFERENCE, "')") + return(get_area_conversion_factor(AREA_UNIT_PREFERENCE)) } } diff --git a/r_app/translations/translations.xlsx b/r_app/translations/translations.xlsx index 5e661811f0199eefe960289bbd3a8a2c6e784a28..0306e69f8686095a0b146d7d6df79be022201eaa 100644 GIT binary patch delta 8716 zcmZX41yodD^zJZp2_hwpG}6rw3JN%McO%{142{w-BOTJx3?0%Xjevjy0|L@5C5*sB ze{22!>%Di^x@+C_?eFZf&%S%@eeT(}4GBy|0_*OuAmqJx%8BRzfC?@EKm-5){G51w zT|FEvU0oe{{G6Stbevs_#E9=gW^S|33zX6tL4DEn$(9JCkM5!~rQO#EY6v2*q4YsCX3J8|JPa!I<@+dVoNFGL-g ztiDk(p(e3?zV z`vjF@{p@@Qb4iOnC+25X&S!P9eyQh2kmL{D!@r2kzuZ_0Vazx*k%XQ}zuP)&`+oCn zGd1mJ0kkH1VSlPsGlkQ;dO(&{`~x%51scKkXd2cY&XuZRLT-z1bk~>aqP3SET^oF% zix7np=FdZ^%KF?g4ibCaoR#E(nta4?5te%q*pn}`$;27x+>^jR%tv)xjnU}^(zsh@ z82V!Ajg}Xm7dkUmBRn=S7T%(@DQ*R!T~4tKH!W`V{dtr3TZke=_M@2A9e*zy?eWoV zCnQvwDB$r5xTY879fTg%I?PKu@Nmd!*B9LU$PdGLb+uwZXx&nvbuLQ73s^$aUxL@T!jm*I`3?j#rcwkcD%{ zRakFDvy^)MV{GS71&p#lA z(Y}oc^>u5e5U)9LkVE&2_TtKDC_SxrQX;52Ibh`}PlH5Fc}?xh8z;-NAa@;m5=rYK7keR;}!I)iRtefKNk zo10y88E)jG(RURLhhpl^?b%2*o{smX?86k=d`ear4$H>C8Y zk$G)bteM<&v!0>~EgNDW2t+7X{(d~Q8>Ng)iwI95hA#cyp4 z4}5t!B<7`JEJT+FopOE;2YSn;Nd9E;KAXL|bvZ?xzgnVoNk7$TG!}|+zM7`Xb4kDd zjpEIp9~5Qu<70P2 zSBllZGhY7t%9b(ZRS7*r%TgD0J7oBD)Xy~6(#I9jfT-F#9E1km(|bHAtL*Z3K19zk zVS#WJ>zf6!GzeZ`=UP+PHn-pGr(dr5T(9o@S!FwWxJsR>NWHs+7orc+;hgiIls>DH!nS<*xI&SMby|B>fhdyI-vNZS9zcZoGM;ji7MBF zm)4XWx817D(o6YV6VOa9$i#L_LrMFLe{smvv}DxCbWk4Q#*S55X*8?S*|pEeg-(wA zcT_=*>J>=GN+;<8sj?bulyT|o5+FA3)A+77q4sKjR6CU$)xhCY;mS94c%h1{GM(sO z%rjs#h^&Hf$?JmCj2A~Eecyl8u8jXGy}&HozkW|NuWu5r!lAzpac{p;de2qZqiA0X zjqgkIlInr?f2XU{9Hm(tTq)MeN`dJsk>Ph2^0-x9;BN zPKa&vlGdrV+KeD(vDW$Ez;7d;k|^fCvzX_Kp`jpfWVG~wRlFL# zlFh(*$z`3lhffLD83yZ1YBh_N}+xuvi_`()J(WC^!oIMRMa z{{$#&v-PVx`x+|ER>Ha>ra}pyjbu?3PHGo{wTXyyvCq^}ph^aNvKKedetE z3qmmcJtTs8B8~jR4`m0I`S*ggm(tQ^ zLOfX92(et%b>{W===zR>4Q@T)_3}05H#<(d{wtL3Wv^&Fp@19*I+khBX?#C8dHVZ| zXAhSSrbV7C$(Cx2K39ZW6gjikdsuxWg(>EG1|(hLU7MP2q*T#sl47ssMmrL&iFa0M zHUpBgpKtdF61kqq09v>4KD2qm({RGv zxNmgLYYZhy#WO6ApVO{qu)PRvTHoa}tTb5U|M3^sWqCX`Pdv;;6WZh?2@wQuB*m~ zS(m)Zai^BAl5$ihdocW5Qt{-xFJgW17J>nNV#mlv&TT_=&v{j0e-}3_XV#TF8_$>zqMHoc&1BgQf|r2*=x zLIy$-2$PYa@yz!L^zrLqlIq~6143_{op!~a*`~j*Yv{+B?+X4d^hVXm{_$hBI-D~@ zD)?eNu=7pCd0$%B>*Aae8uQ#U%crLZQDaP8S-kL`bS}XD@~&f@fMMVe6Rptt(oEr1 z!)r!zm8-@MEV0czLo{8>R4gJL?1bOVjW~RJB{n<+GjrIS!HdIEMdM^i36JMKr>GG$}ou*o!XGzo=KX5G)e8sXlYBn4${80IXA zzs!4N6w|>aT9X!1JqSWda#czCnum{VPGQCJ8j( z4n_0CtH5xa$W|Z&QS6cth`rerNEWc|#^W#hsMwR`+Va2;4huc|7oMuy2YQS?Bpso` zqMz`p2T91$A8EGJq?t^HkXf>+;pv|-Za%<;RDWS(PDt+hh6tvmsdd^rxr%F{Dn{~h z+Q0fe{?o5yfGb8OVYeOH+jjiMjP75C>YwCr=n1l^eeIY+{Q{d{-V&TBOjGP7{w}4# zEGlD4;=?;Rmp?QQi&paZrwJyY3mW)5W%`V~?bNg9hcVa4-*hzgkOdYRUqZ$c8ZmM} z3OCb;^x?^fzwU<8Pbx8NLDwIAv<$E`!xJU??{%^X5}$!jQ~C0#UmlE|WKnF&nwkSj zHn{8iK-tH`;!a76-@nNpFo&ugBXUhiB<9s7%=Yifr3m6o&g^QQnq6b6IqEye4 z+q;BacT{QtYw&K!|R{gz9eixfAOXYz>6r8X~t zPLLRdXlSBmNkt;W!>f9GzArzohV~|7J6|?t z7+siZgD!4PvREduOh_~tIs&CGI((IV;I^MMf;4Nvm!fo7FUD6R+HOw*76Yg1?JI9u zFTZo0OZju{t9S<)q9A|n+!{_go1Z6UTbj{tBJNp@?)#BbFOXex6AM0U*4tq*RvnkV zCcY3uOc|r2p9OZd7ht+Dr~NH;uagFUJcJfgoXYA>9Z?~la>OeQ+S|Ao9nXZMfR1X# z%!SQsh9u-5nCVIm`?RHFZ>}p)0=Cd9K{~N;wp&_7#O8LYe^_~BbUb(UXmf9N7VZ)| zEVX>~d_!)xQ{==9<9_{W>>NUGva9|<1qqp7a>2@rP#WJ@GQ%4C$d0YHf$AP5+>l|) ze=Zz~&^fnagugu@v$s7Z0x4nJR&->2YeQ7s9>?Wu)1Pe)XG`LpQA<-?iC)DKhuGLt z%3U!+%)fWvArDK#AM5&V3VU>19)IWDUM4znLE)a;vX~iebd+svSeZa{empT|!+=m! zFOO+X5Y&by*>e00>An;cF1bI-?bCCkB%qr(K3<`}dHSdARK&w5=L+e1F~(aTfbgy} zCTZmr=GV=3ZF5q}NO4;;fU(zFcVbhVK#tYROE3j1{0?83Rj9NWp)^1n8Uki~byuUF zOW$I;c_KTTx*H=?c8vPAjg+MG`O-lnqQRz!Weftd$zOF{l$!!an&f! z_B&5Rc~7oGDB_-OSkRCBrVCeMJHO4QFmd-?n&YKY_IlSOMaJn|NA~;XUNv$>tNF^D zpk<*iRC;x=rB^Us^Lj`C8kKEz^lNwtl1Ds~?kM#yui}Tg z%eq&Ix*iqc7p933 zj^G%=(#0|g7Ufg@tAtIc0H}p%$}sdG_tfOJt}gk&x(MbZ-h^m3^EFMNvmq@@y9R_@ zRW!rc*~*^{&&@Zejt(kEy|>@WDh7Uytm>xGsA0n$#d5IStLsxvcn>~)JDQ%yEHv>< zan>}Eec?SJ_I|M`K_!6-wwCtC<`H)4kIm%QmZ~=;Ro-_6g?wPlLR_TNT62K#k34Tp zJ{vq0*%>tvFeapMm))Cjs#$WM*gmEw;`3S<6!Qr_HghrVGP5nFIp+>=3Jh|ZCa&iW0Gm|f6#%Yd1TD=@^_FbC57y+l zLe61!B=|D{=Q%hn&CK-r~-XrWeVEfM)@w9JhMl>YJ<=}O?G(| zR{MvWgC^bS1@h2M=%{{lksTZ)G&W^Puab00YTmO_$fTS$8M^4G8N_~DIH$&HMCdebY3?KUSwg>;n( zzfr5IKl`4Fx1WAT_q~rVS)1z8W{gT#CDA70R& zK&yA=ZzP$9`kec1_}5Fwmu8nse+ByN82>=0(V*GqhrXV!aGY(_XHY8}e&rSk^@WcS znYAv38tT_{ZoJSY&+(W0OE{?cc>Z#bW$V)a?0VNvXAZ@b*;?zPb5$a`2SQn1tgWew zUf(8e-EzJ)6{{t@J=>Xr(3jqa2+VYGX0=zhAmh2F)TGtQ=%*kX{-y+JrDxq=;*-xx zgKFl0iKU-iJ;4|e+^#}fkRrGnUuoaBs=SH`#4TGwVF%ZTUL6|L-3R6MpPabsC!Bvy zFL}NISG!GOA;|1~6}jhCHk(45K8Gu@I;`lvA)ZLKX?mUO8HeNC8tFS@myX8eLm1yC zd~_PsepmJS#Wds#z5pVvqh7sj{pa~F^G>Z$HzBScZ*`# z9@#Gr+s*dV^GEmA;6>b|)5bdWRU&5f_1x@A_B)vi9j0><0y6g!C`D({(fU74kl=C= z<3b{4~}$`f%DtwCFai#f~GLB(~ZM@+ax(nCMq^u0V`wl>Fc z5#YOP&u@gprlY`pMfn>8sF!R6Fk%%fo?bQdFYDm_ptgKw$(#3alRbbxZ_P16s9tWd z9fpqGGOzWj!F!^%+2vscK-Nrmr2DgRJ)L@<6$NI=-+1rSb||@ADy)E|8Sc}`g-X<- zV10ZJr}W68*53qV#3%7Z8_N1`bZR}_-+=Yc$gUty6c~S&Cyo;PO3#v6k!IAw6Ka{O zu44mZi!y)@11{o=qCF}W7Ynzb=R4IsUOGcUl*r9G?EO?UN!1uky~PyRBs=lP*O^qy zQCTnnx8;|?;7I}o-?`ZG$Vm>_P*L2TA8Ml(L~-Oj8};rn)$B?;Kc7if!2wL2z<3t4 z9}xbO)_&b(7O%9&E6jZ$8grytTXWI{c|0{R}omz z%~>n5IzyRn-r$*n+hTq(#%v<%w=il^zS_EVAr`6VX@x=nB{e42m|JT9cq@%XUcc@D zE9L!fXtPmF1mzEyJvN*JT)m^jmzDw*fb>|ra~o~6=HqBkmf&l|$|H;7uyJ6=%~3DCEPzJ5AO4hWRP!d`ih z(q?M*#g0nKTT2U|_k8K3NrN^#6Xa3hu#mH+E2?KnLkoqwxcl`H0H_0uFUGnDRhjG1K!9 zr<@?38ss?@I|@6m+ae)~!i~7a>F}(fRS}^Tb+Pl4dS0qElhqG}iQ$(S2hJ3rJFz=< zkcsxPyxp%BL`nUBiHriY zF-==A9o*F;U-{hl;u!I&c)~N>^%2(u5RJo&O82CyPvf@3HIyCMXKr|P7WpV$602Di~Bw(ZxU3Z)y2jq9A=uL{K_zl+$e-Gx}s4QUw+C7>YDTJf%fI*NwW9wZiz zthxuAUP-cf=ama%bUWPq zsT<}nQM|cz<8a$PZPAw9W=A$%0v3|KMlmJ-RgQc(dxp) zcF0$@M3OWzF9WB`71^zg} zBD_yZ{_CMzUXH8XWePYqt|L3a$gkf}s4|5Br;=jR{hh)18Ec~jBnrDdRnOFU*9zE0q1Nnu7i%;~G=kH@!Ogww518Iol3yOeCpRU;XqlUR~Ry znbJ^~sXeAqCOu(RcEPMgf*hEbO^2jT;$uhF&XZ8|@pIz3xrl9?Ie)im zgg^pxBB;4bWz4HjDtPTJWt6+*c*YxP6aRwNNTUPs+rv40Ip7(6I8Hw2pIXW53{Y2J~Wyc zak!e@Y;>spS$Rd2XiA*fgy!wjpfGr$=F?|IIG7*x2cr-+J5Bo0o}~H;WM$SR46hXG z>{86kN|f2E1Wued35Rg&!q_;n`Em<#Jwd&T^|Wss6uJt6TFc{$Cn59`oU_^ z@Z$OEKV720-$;Jm$SO7MI;^!z?LE?a7e5b)H!;|Beoe*HKXa;vuZRXoHr-#BS_t76 z*~@I93Gtl&kmXC;P(<}O3gE@)=z(B{|6h!&Q1UT delta 8765 zcmZ8{1yoegxA)MU0;0gsU4wKe4blxmmo!K>GJr^jGIZx4Fd#KFND0z1C>;_?3rI5# zeCYqZ_kHiXYu&ZZ`u+CW`<#1r+;jI$TfuB!!K}Z)uJ(OV(y#p~u)nF)bfZ z&PfVDOto-th$wD`Ysl^_>K;LIg)hn@Bh84<6h@7lb6I$wzxjk?{_Zgw(GV}6?MN03 z>>cH+N6u1%Y~~*m=hbw3<;J|VJ#UkUGx5B?5ypa8e!N~jNclZ6s!qwNc63t6#P;?p zJSoYD7bH~4o*`JPXIayD9AYY*vd!MsEnnCISwB8qPvxUM_@%)lb#16DJc(`onF!?Z zQoiUamG41`3G;mZ(03V)x{(oKBgy5mjazBG>rnICKHA3H`Nx>92nzz{prhuWhYbri z-+n}8ZrN&Bl6){O@~(Jdv}YFohW>|PT3GR=&B*LlJ%-hZBeiHVcGb$|p`7@K3%{Pt zdPu}FnTA?WGJCu!f26DQ4n)0XHoyYJ)F8%jCz}4?Jwq{h8y4%Ng3QA}Fpb&3SR^oZ zMi5gbU)`)BWbwpG^C+_Y$&hG{d~F(cQRJP$h&O)#0`#<$^I7eiH72e}!!XUb5t_>j zv_}}=8JRTYs#;wi%HIbozHNBwEzR@r9!SliTb#6yqyLbi2r5SR#pQ|!sPrp>MBE#* ztSHHjRi)*}>+u%5ZE1-CG|0Z~-6w&~(KER~=5EeX^R8XC&xCZjj5c9S0jVWF#T7YG zb@ayFhmmOr2`@!{!~|=tZwqioJ0k3Oufm+JWs^7MgZVYpD zB*~M2p5TDoy$KDD+}dz~?&J`kA!lA9KawiJ)Smv*ULt6Yzi3+<9Rvid(x7c3!( zJBf5k_=;~h=$&-Dv{FV``Z2jyW2UjYhHY2{>Jvw$r{ia>Yb#Y@pO4HZ=j(_4FZ0O)!s zRZ}~m0w!t9nbJECLIA*B8UUaLKs~Fep+FKEh}sBg8lN#G4MQ;w%}T(f=_XgaP|FYX z(LXNJM^?xHcqr*WUHhkiv>e8xZ`8`a-)3#CkgWwQEo8MYimlv+7p}6Mns!z|v}{_K(gD6ixrO=Z=;ZmE2ky^V zb~5T{-vlfriQTR{_|mR+9nQ|REs3JH=LTZ^TcyDdi^P^DuFyMv&z5REWuF~&_FOk3 zj?nzyCPFT1Tkjk$?)DCah9Q@Q%}a0;gDMj+xFe7dRc=@)hFnap<8K4Tj=t< z;K8cvY^%=YDdWsiNY3~}vcO^7vwav^DsIimKJdy7C1QLS;6}uYQkqkc^b)}2()e@fWP5I6^U6nd z>Nj0Mj$O{;@22y8P(bIJHrvxBND<%RoTS-ek!R-K5LNumQ$=$@RYen7r{NL*%o}$X zNY7KBVCWtjME2{aUCY#OU~&6Xuem^LrrAw?M|Y5g z`d$2s#`<6wTx7~^8U;_|Q#RJfw4Cpc35=gMu1T7HHK(~G_WSNN!Z^zYB8sywW5!vk zRH$o{QO+^O?XH#u4UJn0|F*wBOVbmbr%mk*zoKON7RFv8W_T%@O><;Yk3OEmb0yM= z)UU_@RL$}dAl735u(3-l#LDrFoX5IG@j`-GM0CPZO*BUy&*7DjD~KfQMutPzJZR6@ zq6|=FX?TBs1x%o*%^-t=E=1NFd@w;j*HS=L3on8HWTbjWKGPYy_(^+6e&ZX#NP$Xm zSscQ*Jsvo3tr~D8x>*L>6XIuR3OC9rVCf)g#qkqOtmwJ+6t3 z2jvYw)W#A!22Y$ZdB<1Z6`3&?O%bXEw(|)3>R+OXXdo0_=VI5W_*u>i3?VQ`tOC$! ztyOaD_gZiyL()V{%nYN%t1{+ataj1Vi=VS>M^-63X7{LT#?LmQ^~xl#1JHeHjGQ%f zrcFV%Q8hdN@Ry2jismG>@t~UO$8a(&j%ng5@vW!~w(l<>^&Zv2tQV6}FZQRQjr|WE zVOp);FY|_NynwU;jJr_9LzOJ2!GACdf9g}|#$Avm9_0`gDZ z$d+LKszBRppM7d=x}6r*5L1>efW>WKRSErz+j)R9F2nu`B<3y)oDYlnmDDon_-gEt%;_@!JTc%@gjI3PpEC&<4z$NhA z8|;eBcg=p2u*+#H4v46LIZ2vpX_m9!5(4Ku8z9@_1DNVg33u(=vY`n57+E6!(-!S89?k zeR3h4NHds%P=Z(DU6ceFS3!!tlu3x|K)xv#!KnoCyjNE3v=o^Kp)z*?ZUt|U0LDCD zaEO`&h3felxaVC+Pmkh6=H|Gpcv3V@mT&w%Ge$df+K(4@}II1_}83yE0Yy#LDCZ{nD-4 zkXdo*t?5|qgjRg4RxvR(T=HL-@c(?Ktmp3P!R0M;u#u4}By11q9L=EVE@CHaaPFvF z5$(+lbtM7Et33`gd?32fiDa?ygnWE_oS+(3KWm5ZUv12sRo6sy0}~c8oJw zHo%~ph_~!j=exakB6$!Kem*4^#p$Xt4%pkVTZNCq<$fouXu;upl813YpV0L68{;9# z;PGJ4cP>ex`WA}P;Kg9ot_@qje~7JD_i7OL#Jbi6%I4)j5yRsEp}b#Y3%$+ z!FE3w9%H$EG97vB?fj7^)oiKles#uDy*ww{?d351V%7_e!(bpKN z+?E0v+QG5s06JM4u*6AN+kH+sd-6|NZd9n6Q2%`(E6?s*41Bj>yRMJfoTS~*leY@} z8PRTH6xWRv7N<<#3-Zf{VHp_LwA$lR6hjd^Vj64mecx=8IJN z@V^DnKp7`Bn0`VxsO^rH0UNr6t#S9;V}?D>DNKHurNS9bI0uU6T1m9;YJ7?G%6i=6LYj=ypQMN6J%j@M8$t5$?%#V)EqYK#DiM5B$ z?bqasm}S?cHLHU|0GiP1?o=2HLnf4 zLoE+|bS~B*|6WB&g}=z72XO%aRo4H03xz!Sdw$H%HTNu&pc-B?Sv-me^6a(wtgNN_ zNfU8ek92ktoJyA<$z^*H*<4_G7`SGSv@05*{lu!pnYDJZ{(bBkmAxDu^z@+HX`$KZ zl?BS3=B7|4^BgAitQIXB4rw#JvRa}Z)@t~W>T?Kz*50lE+Aex&z4z0t6WVy%3mIsz zM0M6l*G3id2maP=Zn|Mhm22_$<&wIZm%%x#diM$D={Jw3M+@sKd7pCoX`s4&NVR*# zxMF<&ci7wEUmNA`P^GW0mMr~U*SZ---If^=e(qoU8O+8$XbXa07+ux|?uEwsobT>} z#XoHVf1xr?O3$q5t{(PT{|V2zMM0)BNR1Z0dUPK;ewS<;l<7e`)DtK89E6pt&n<^7Ddq)R-33Of237zB~P6-}5I3s%jfPu)T5nGn}G- zn2)4v=h+q+A8P7j=ZYbyf3L1f<>lP%Q}ZABsJ;~s*~nY?mWEr`<=3z<_&($VVOZ3D z^!L(?P}Az1b?W2ad~cSI=E7mHv$bx6AR^M&eP!3TD=0n3okCL$0U%Ssnz!lmaq|ls^#g%kP-Cm!uhHW}vc!x!d!KMZZo2V}&q~{S*f) z^n*;R<*TcGg?V$7h6ZhX4?Co*2eEZ(+WZ`xR-KvJcnIX7o+y7v8?{Q*)f1c(+oYU! z>?ibbl?ZZ*`Qa$AAzKJ|@h<6zxy?^u?PKsXbzrL&dZSM%7QAFog2XkCyq94)D{Sn> zvv2FB!9BYc=piYzO&=A{%;6i<%uce+%;S3u2CS0^DL||lS!v0HAS!Hc$D>bt)cQqS zCG+w+|efn^q#uOg55? zuX#>>X}~jVxQ4%H1g4pbB^;9a8XfJb0QM?3p(vQqr64)qk}rV?lteXb zd0<;&zMAJjCA!-JkblIPxuZfM(&5ccejP+92gYT{PG#c64FW1S5LlXU{9v{;iJNq3 z-Dv4}s0-J~0{6|wWsJ9tmxpOL@xWe0ZN)v&aCE32^Lu`0rYoi^Z4gWC=&Gq1!fDWR##Z7& zqnCT(qtveLE&=8|w*xykOJ$d0BD7Q%00_kAhrUQbhuTvz>NjdUd+{TxJ+#y(-$-PFU? zNEd(fIJp{!I`en=sF~roL_GM#sG)+t?fs+;kYdAM8;=r})ddaZl5YYM2>mNcrd-;g z_O84S3Y7SsDD#zcfbLZ{Cwvw{#DTZiAaBCU4x z_0X12J4{~kt)zdI@!BM1wLt70IbS8}Ng~QQFK56^`ATKVu?6O%-{;jf(-8u*nSPzZ z_f`xg@LL=RI13W>OsSKRi)=LI=UikChj+B~#z1TA{n)MF*1Ajp}jF~3-B!k5P?3vMQo3Wp<(%g`O##xL?zf(ofBVTzC zHE9&TC~ly>l!FJh<~FGB41=&wM8n3N3dec6&(|v}))|ImFYE5~fgQ!#$%8EYuxB1$ z|8_?I-nlfC8L>|MStZr7YV~qGcb)jjOW)wk)sIJe26wd~OqPMia(~b>q~Q`Q>!w~3 zuUFS*{V#7ZCM9*#*904-i@FTF(E zpLyIUrfH+QHhri4YE5RH&i@ql`Tkrt_h2|zIwdT}G4r{p3P=PCt9``@YV7 z`|mpU7>WRrL+CK49cyqp_S9W;ml#(`LywfjJd7CUp5^g}=>W$<3>b=LZ;w7cT^{V23;5)1aIjxyfwF$Ff9G|H z4s8A4f7-sHwh&CdK3|8noGUOs4L1L^RpRV)t_H+dfXq&m?0ld5>9Sh**wj6wav&VY zM>G9L*HHD2-u9ctgR<3^HTpmOum9FH`oDDz>-Mlt3X$a;`#aYww)FZ>SO@q10D*#x zF{&Dc#5?a+Y`;17gXe8{lWm6EmvVKt;@%Q{XZ$DwyX!nDS2PmyH}~v(q*PcWnKz*G zbEFHx75n_>)d#|X!LE0L8D_mghQ}W1V2MG(A)7yRuay-+qLveT?Tc_Uq){ULvoOUf z8^N|Xixx6ZM6)z+XBZU=B9r`R-3UQCN<^Rs?aa0P_Pq#d0nfm!C zm1+NOuc#(GCU!Ema!=NpR_4s~w&6s~$6V&9(|Uid=*4;SeuO%{y~#+^V%Ug-xx#FJ znUlP1dgH1h{?*68_g^nw!+BS4rsiB{4rjy3CzQwq>~Sxwu9pMO*QbR0scQNj7}s1v zqNfr$i{6ak+{vnfb5<#ytub8N(ZrU-ZTvW6{nW`BW(9mSl?_(1u9T8ZMHCHn?Id(A z!NSj(wr#&0R_Rhl*D}q{Z4<^+yB~e)uS3%q|N4=#Hfe1mG#beFtpR7~34=ZDm_HwF zz2p}<4$VbXttgpu#jwHwjjWXah zmlU*Q7HXvxQYtj!GeRmc+4f)3N;Mwl+Q7BA7X-Lfb~jkl-w{j+Nu^+GRqeg1gB0q6 zX4{bHXuT&ta#9lFeuZs~4VlVWJa|}V1NW{7B`g!cAs_t=-=;;R8&nLDiZa4gTasWk zgs{GJff%et0m;bW_2LnSmn;>$ggM2AZyxWY<)$-xEz7-9o8)KW*ef^GoTD1T&r9!( z)`?k8u^aiJS2tt`X${XpT#&n$hJi8tmAY|zPwi&M{jcxwrG@myU zq7h|ISm!FoAV!DGje`eQM2rj?gfeKV+Iea`1khl zkB;SwXql!Oe6G?Z&lV%Zr+N2ofUAwxHh)07f9UP4SJ~n0I*8dTrfXm@6H>1I`_tP| z{4!#j?D2a@**V>Gi^0t7b?<$JL$rA8`K zb4XF&T7OX8d-&x|c@k8B$*8Xxr?{IXYM}R?wKrFnP&17vPcm>|0)4Y0CfbC@MG^fM z=&m*fh;&)$ht!u4W4D~O?d9Rb`by~2$!Eyp&>1Lv=}@)_8AbG8g<^4X6V=l`l9VBE z8q5O*)`B<=*c>A3#G+aRxub$LrF{<=47geQ_zQx+vpI~0fPLl3-uxvx?+Lm*QoL4? zu1F}7+by~N9VOH&raPImpf|H|=x^$Gi*_-E&ERI0SZWG+ql*@;qDwat~a=u_3x}gUgLxIUs^e)yl zwkq2CL(>hsB`Xsm`4sbwQ|(Y=ZZ{9TWDj8ir}-rbPRbxdDVaCgW)&}=i}Vk*duEfk9=94y z?rCo8gXBpWd@dc@v&LnUUHGW}64_?dyP7Hd+6&3!@Dzsq6jqK6J0yz98al@$ed*Ho z=6)*s?~c%fI%hA`kBHQMuTc@Bm#do~eui!4yZ6P23 zw7i}JvstJLI=&TH>i8w_u{WbD0i_l&#d&0y|KhajhHe!YxnU$uZGpZ#s)C~pR{XT| z{ExdbLbyw>>vxY%*V$vb|A0j-_cP9$LJT7Mj4qGorloe))MzKxfnh<++rMQ0(|4$5 zN8m#C6T=M^?;D4i!~g)ve;S9re!=bzzJHH>2PPvgmq6s*4|inAyLip)+FO{p<2Lc} z-vvP>&Xu;81F3ky4S6e))`sd)KmAZ*oP5;6p! z(nei&XdV^B>c=WJt9?8lfk_%HNC!glY;3hhEosHKV_DyoLhd)JsMag+mhpsQw)Ir=7=NkstGxZnshkSI{ye2f zXfGQF{It5`#F8AQ<(Z?wIGLud=B8T`fK+%A`5tN$!1}FVe1&)Wz|kSTU9~m(HRoV> zCRGA&GG94qy4%w<{KbhjruX@n+I;Vh~E zdQ;5ZCc`m$Vn5ms7ozThB}geqcEGBqqo>G1tey<`dn=veh}kYBPv_ibUj}==F-@Fn z%R)z`k?5bg48&VBANQjLF10qI3M@%dru-6jq9vW^xSdUO<*wE+9SOcdjAK)OiBy?1 zyYS9&&lutb%e{@|KKUrQq35ZTz@-jUGv8S_wQR6*9qlFDao7!w!jP{T<4@mG?=x7%FW&PBNc{zt-gChmLmAb2HBmYrA*ziWXU~Gu>CtYoFSr z)&e{1)OZ^=1N_5sHp2fDX1a5hZo*R!j2CEFKB2D`3oh2@MGbf6=5Xt{j+LY*^#2_C z+`jB+8)x;T9Oqo!NgrAjDHQAq@)NS=`|U!G^KcC>>+?ek<|S#te4FxhXGaMOk#pvG z#kjFmJUlZ`zq6%HfHPUUPMNnCFg7x)2)9K1F*IKMQPGeUFm~EbCf9l>SmrxW_O$xZ zmF>Kx8N+puV@@%>cEj+QiJ~AI!MHW{y6lVoqdVOf{}pM&g6Y0E_#i1#lO~%o4h)E% z#sG-w%^J>~;z*NaZU9|VGsF9Xk$4;b4fE%}TIxAS)S8z*ihHznl;8)eO|~;`Tck}~ zHWK#qg7s?X-|0nT>cW32c=()|3G_cAgqJ1ER2#o7!vc#4)_r&aGccwL5<}o_V$b(N z)3tIYr)MP#?A#Y<*buBox>go*D*BM^500-z9i%%j=vAkwG`!LcRK1Ca&*hX<*_^QS zA-UC9UNI$@1f!;!D%f{H474tH^gHRol~liqUsc*K`-(YrOSD}Y6A!^$RJ4$vj;{&A zt)KPy87Y=5S_;W251Eb<>zd?kkZ5A*u3+ofWvvUBSlTJqn%}(ofOFibs7I)0T4RIL zfB63Equaydt-9x_c)%-X7wfHZcp56!Co(~T1q#ZInJ#B*N?(SJT$Pb{39YFdGHwYC zeH@(M+yDWyqRVeJ(}Xt&bp6@|KdPRP{OhA`SD0fF{LhaNP>m{WjAm$V6+iv|?ycS@ zx&OH?i2rq3ud1l&|M#b_|B<*K