From f879cc7d3d9fb6a7b3bc8884e607ae6244b633dc Mon Sep 17 00:00:00 2001 From: Jeremy Kidwell Date: Sat, 5 Mar 2022 08:44:38 +0000 Subject: [PATCH] final fixes to conversion over to sf() and tidyverse --- mapping_draft-hpc_optimised.Rmd | 138 +++++++++++++++++++++----------- sbatch_wild.sh | 18 ----- wilderness_layers.qgz | Bin 11202 -> 0 bytes 3 files changed, 93 insertions(+), 63 deletions(-) delete mode 100644 sbatch_wild.sh delete mode 100644 wilderness_layers.qgz diff --git a/mapping_draft-hpc_optimised.Rmd b/mapping_draft-hpc_optimised.Rmd index de88e33..e468eca 100644 --- a/mapping_draft-hpc_optimised.Rmd +++ b/mapping_draft-hpc_optimised.Rmd @@ -386,7 +386,7 @@ tm_shape(admin_lev2) + -```{r 02_admin_ecs_normed_choropleth, fig.width=4, fig.show="hold", fig.cap="Figure 3"} +```{r 02_admin_ecs_normed_choropleth, fig.width=4, fig.show="hold", fig.cap="Figure 2"} # Plot out first figure with normalised data: @@ -403,7 +403,7 @@ tm_shape(admin_lev2) + just = c("left", "bottom"), align = "left") + tm_layout(asp = NA, - title = "Figure 3a", + title = "Figure 2a", frame = FALSE, title.size = .7, legend.title.size = .7, @@ -414,7 +414,7 @@ tm_shape(admin_lev2) + tm_shape(admin_lev2) + tm_fill(col = "ecs_count_popnorm_scaled", palette = "Oranges", n = 5, - title = "Concentration of ECS groups, data normalised by places of worship") + + title = "Concentration of ECS groups, data normalised by population density") + tm_borders(alpha=.5, lwd=0.1) + tm_shape(admin_lev1_simplified) + tm_borders(lwd=0.6) + @@ -427,14 +427,14 @@ tm_shape(admin_lev2) + align = "left") + tm_layout(asp = NA, frame = FALSE, - title = "Figure 3b", + title = "Figure 2b", title.size = .7, legend.title.size = .7, inner.margins = c(0.1, 0.1, 0.05, 0.05) ) ``` -Given the way population and places of worship are unevenly distributed across Scotland it is important to represent data in terms of relative distribution. For this study, we attempted to "normalise" our data in two different ways, (1) as shown by Figure 2 above, by taking population figures from the 2011 census (see data sheet in Appendix A) and (2) by adjusting relative to the number of places of worship in each council region.[^15914204] The latter of these two can yield particularly unexpected results. Thus, of the `r length(pow_pointX)` "places of worship" in Scotland, the highest concentration is actually the `r as.character(admin_lev1$NAME_2[which.max(admin_lev1$pow_count)])` region, with `r max(admin_lev1$pow_count)`, second is `r max( admin_lev1$pow_count[admin_lev1$pow_count!=max(admin_lev1$pow_count)] )` (`r as.character(admin_lev1$NAME_2[which.max( admin_lev1$pow_count[admin_lev1$pow_count!=max(admin_lev1$pow_count)])] )`). Rank of Council Areas by population and number of places of worship is also included in Appendix A. +Given the way population and places of worship are unevenly distributed across Scotland it is important to represent data in terms of relative distribution. For this study, we attempted to "normalise" our data in two different ways, (1) as shown in Figure 2a above, by adjusting relative to the number of places of worship in each council region and (2) as shown by Figure 2b above, by taking population figures from the 2011 census (see data sheet in Appendix A).[^15914204] The latter of these two can yield particularly unexpected results. Thus, of the `r length(pow_pointX)` "places of worship" in Scotland, the highest concentration is actually the `r as.character(admin_lev1$NAME_2[which.max(admin_lev1$pow_count)])` region, with `r max(admin_lev1$pow_count)`, second is `r max( admin_lev1$pow_count[admin_lev1$pow_count!=max(admin_lev1$pow_count)] )` (`r as.character(admin_lev1$NAME_2[which.max( admin_lev1$pow_count[admin_lev1$pow_count!=max(admin_lev1$pow_count)])] )`). Rank of Council Areas by population and number of places of worship is also included in Appendix A. ```{r create_admin_proportions} # Calculate factors by which ECS representation exceeds rep by population and total pow counts @@ -450,7 +450,7 @@ Turning to the total of `r length(pow_pointX)` "places of worship" in Scotland, Whereas our initial measurements indicated a prominent lead for Edinburgh, by normalising our data in this way we can highlight the stronger-than-expected presence of several others that might otherwise escape notice because they lie in a region with significantly lower population or numerically less places of worship. Taking the PointX data on "places of worship" in Scotland, we find a less dramatic picture, but also a slightly different one. The positive outliers include East Renfrewshire (3.4x) Edinburgh (2.9x), Stirling (2.2), West Lothian (1.9x) and Aberdeen (1.5x). Again, negative outliers are far less dramatic, with only Midlothian possessing a ratio of more than 100% negative difference from the number of "places of worship" at 1.5x *fewer*. -```{r 03_admin_barplot, fig.width=4, fig.cap="Figure 4"} +```{r 03_admin_barplot, fig.width=4, fig.cap="Figure 3"} # comvert admin back to dataframe for analysis admin.df <- data.frame(admin_lev1) @@ -465,8 +465,8 @@ admin_barplot <- geom_bar(position="dodge", stat="identity") + coord_flip() + labs( - title = "Figure 4", - subtitle ="Comparison of Groups by Admin Region", + title = "Figure 3", + subtitle ="Comparison of Groups by Administrative Region", fill = "Groups") + theme( legend.position="bottom", @@ -476,13 +476,13 @@ admin_barplot <- # TODO, need better export method now that tmap() and graphic device are fixed pdf("figures/03_admin_barplot.pdf") print(admin_barplot) -# dev.off() +dev.off() ``` ## Other environmental groups shown by concentration in administrative regions (NUTS) -```{r 04_choropleth_others, fig.width=4, fig.show="hold", fig.cap="Figure 5"} +```{r 04_choropleth_others, fig.width=4, fig.show="hold", fig.cap="Figure 4"} # TODO: consider switching to two-dimensional kernel densities instead of dots as shown here: https://github.com/mtennekes/tmap/tree/master/demo/LondonCrimes @@ -503,7 +503,7 @@ admin_lev2_scotland_ecs_plot <- align = "left") + tm_layout(asp = NA, frame = FALSE, - title = "Figure 5a", + title = "Figure 4a", title.size = .7, legend.title.size = .7, inner.margins = c(0.1, 0.1, 0.05, 0.05), @@ -530,7 +530,7 @@ pdf("figures/03_admin_lev2_scotland_ecs_plot.pdf") print(admin_lev2_centralbelt_ecs_plot, vp = vp_admin_lev2_centralbelt_ecs_plot) # TODO, need better export method now that tmap() and graphic device are fixed -# dev.off() +dev.off() # Second plot, revealing transition towns @@ -550,7 +550,7 @@ tm_shape(admin_lev2) + align = "left") + tm_layout(asp = NA, frame = FALSE, - title = "Figure 5b", + title = "Figure 4b", title.size = .7, legend.title.size = .7, inner.margins = c(0.1, 0.1, 0.05, 0.05), @@ -573,7 +573,7 @@ tm_shape(admin_lev2) + align = "left") + tm_layout(asp = NA, frame = FALSE, - title = "Figure 5c", + title = "Figure 4c", title.size = .7, legend.title.size = .7, inner.margins = c(0.1, 0.1, 0.05, 0.05), @@ -597,7 +597,7 @@ tm_shape(admin_lev2) + align = "left") + tm_layout(asp = NA, frame = FALSE, - title = "Figure 5d", + title = "Figure 4d", title.size = .7, legend.title.size = .7, inner.margins = c(0.1, 0.1, 0.05, 0.05), @@ -607,7 +607,7 @@ tm_shape(admin_lev2) + ## Cartogram Comparisons -```{r 05_cartograms, fig.width=4, fig.show="hold", fig.cap="Figure 6"} +```{r 05_cartograms, fig.width=4, fig.show="hold", fig.cap="Figure 5"} # # TODO: plot as animated chorogram: # # https://www.r-graph-gallery.com/331-basic-cartogram/ @@ -636,7 +636,7 @@ So why provide this kind of data (i.e. at the level of individual churches) when ```{r 06_ecs_denomination_table} # TODO: fix to work with new sf() object ecs_dt <- (st_set_geometry(ecs, NULL)) -knitr::kable(table(ecs_dt$denomination), caption = 'ECS by denomination') +knitr::kable(table(ecs_dt$denomination), caption = 'Figure 5. Table of ECS by denomination') # TODO: Add dataframe with overall church numbers for the UK from "ScottishChurches" dataset by JK ``` @@ -694,7 +694,7 @@ The key question which this analysis seeks to answer is whether ECS, or the othe Of all the groups surveyed in this study, Eco-Congregation Scotland is the most heavily concentrated in large urban areas (33.53%), exceeding by almost 50% the rate for all places of worship (22.96% in large urban areas). Transition is a much more modest 20% and development trusts a bit lower at 15%. It is interesting to note that the rate of ECS concentration in these large urban areas matches the level of overall population distribution (34.5%). On the other end of the scale, Eco-Congregation Scotland is the least concentrated in remote rural areas (with 3.93% on level 7 and 5.44% on level 8 on the urban-rural scale), though again, they correlate roughly to the general population distribution (3.2% and 2.9% respectively). Places of worship outpace both the population of Scotland and the footprint of Eco-Congregation Scotland, with 14.98% in very remote rural areas, but this is exceeded by transition at 16.47% and both by Scottish community development trusts at 32.14%. So while Eco-Congregation Scotland correlates roughly with Scottish population distribution across the urban-rural scale, it has a considerably more urban profile than either of the other two groups surveyed. -```{r 07_ur_barplot, fig.width=4, fig.cap="Figure 7"} +```{r 07_ur_barplot, fig.width=4, fig.cap="Figure 6"} # Need to flatten urbanrural based on all the count columns and generate using ggplot urbanrural_gathered <- gather(data.frame(urbanrural), key="group_type", value="number", ecs_count, transition_count, dtas_count, permaculture_count) @@ -706,13 +706,13 @@ ggplot(urbanrural_gathered, geom_bar(position="dodge", stat="identity") + coord_flip() + labs( - title = "Figure 7", + title = "Figure 6", subtitle="Comparison of Groups by UrbanRural category", fill = "Groups") ``` ## Eco-Congregation Scotland\nconcentrations in Urban Rural 8-fold classifications -```{r 08_urbanrural_ecs_chart_choropleth, message=FALSE, warning=FALSE, fig.width=4, fig.cap="Figure 8"} +```{r 08_urbanrural_ecs_chart_choropleth, message=FALSE, warning=FALSE, fig.width=4, fig.cap="Figure 7"} # TODO: Double check data licenses for tm_credits # Generate code for inset map of central belt @@ -732,7 +732,7 @@ urbanrural_uk_ecs_choropleth_plot <- tm_shape(urbanrural_simplified) + just = c("left", "bottom"), align = "left") + tm_layout(asp = NA, - title = "Figure 8", + title = "Figure 7", frame = FALSE, title.size = .7, legend.title.size = .7, @@ -794,8 +794,8 @@ simd_shapes <- st_read("data/sc_dz_11.shp") simd_indicators <- read.csv("./data/simd2016_withinds.csv") simd_indicators_min <- simd_indicators[c(1,6:17)] # Combine shapes with data +simd_wgs = left_join(simd_shapes, simd_indicators_min, by = c("DataZone" = "Data_Zone")) # Original data is in wgs, so need to convert from wgs to bng to work with analysis across data sets -simd_wgs <- merge(x=simd_shapes, y=simd_indicators_min, by.x = "DataZone", by.y = "Data_Zone") simd <- simd_wgs %>% st_transform(27700) # Create simplified version of geometry for visualisations, with a tolerance of 100m # Workaround to address invalid spherical geometry issue @@ -806,19 +806,21 @@ round(c(object.size(simd_shapes), object.size(simd_simplified)) / 1024) # Augment each dataset with relevant (geolocated) columns from SIMD # assign combined table with SIMD columns to attribute table slot of ecs table +ecs <- st_join(ecs, simd ) -# TODO: need to resolve issue here with simd columns not being added, also issues with geometry -ecs=cbind(ecs,st_intersects(simd,ecs)) # assign combined table with SIMD columns to attribute table slot of ecs table -pow_pointX=cbind(pow_pointX,st_intersects(pow_pointX,simd)) +pow_pointX <- st_join(pow_pointX, simd ) + # assign combined table with SIMD columns to attribute table slot of transition table # TODO: fix bad geometry here causing fails as it lies outside simd polygons -# transition=cbind(transition,st_intersects(transition,simd)) +transition <- st_join(transition, simd ) + # assign combined table with SIMD columns to attribute table slot of permaculture table -permaculture=cbind(permaculture,st_intersects(permaculture,simd)) +permaculture <- st_join(permaculture, simd ) + # assign combined table with SIMD columns to attribute table slot of dtas table # TODO: fix bad geometry here causing fails as it lies outside simd polygons -# dtas=cbind(dtas,st_intersects(dtas,simd)) +dtas <- st_join(dtas, simd ) # Augment SIMD with group counts @@ -829,19 +831,19 @@ simd$permaculture_count <- lengths(st_covers(simd, permaculture)) simd$pow_count <- lengths(st_covers(simd, pow_pointX)) # Generate simplified dataframes with counts for each group -# ecs_simd <- data.frame(ecs[c(1,16,25,30:36)]) -# TODO: adapt code below to include equivalent rows from above now that we are using tidyverse -ecs_simd <- data.frame(select(ecs, row.id, col.id)) - +ecs_simd <- data.frame(select(ecs, name, TotPop2011, Overal_SIMD16_Rank:Housing_Domain_2016_Rank)) colnames(ecs_simd)[1] <- "group_name" ecs_simd$group_type <- "ecs" -transition_simd <- data.frame(transition[c(2,5,14,19:25)]) + +transition_simd <- data.frame(select(transition, Name.x, TotPop2011, Overal_SIMD16_Rank:Housing_Domain_2016_Rank)) colnames(transition_simd)[1] <- "group_name" transition_simd$group_type <- "transition" -dtas_simd <- data.frame(dtas[c(1,5,14,19:25)]) + +dtas_simd <- data.frame(select(dtas, Name.x, TotPop2011, Overal_SIMD16_Rank:Housing_Domain_2016_Rank)) colnames(dtas_simd)[1] <- "group_name" dtas_simd$group_type <- "dtas" -permaculture_simd <- data.frame(permaculture[c(1,11,20,25:31)]) + +permaculture_simd <- data.frame(select(permaculture, Link, TotPop2011, Overal_SIMD16_Rank:Housing_Domain_2016_Rank)) colnames(permaculture_simd)[1] <- "group_name" permaculture_simd$group_type <- "permaculture" @@ -865,7 +867,7 @@ allgroups_gathered <- gather(allgroups_simd, key = "simd_category", value = "ran ### Barplot -```{r 11_simd_barplot, fig.width=4, fig.cap="Figure 9"} +```{r 11_simd_barplot, fig.width=4, fig.cap="Figure 8"} # Run plots # Faceted stacked bar plot @@ -889,14 +891,13 @@ ggplot(data=allgroups_gathered, labs(x = NULL, y = "SIMD Rank (in bins by quantile)", fill = "Groups", - title = "Figure 9", + title = "Figure 8", subtitle="Distribution of Groups across IMD Domains by Rank", caption = paste("Jeremy H. Kidwell :: jeremykidwell.info", "Data: UK Data Service (OGL) & Jeremy H. Kidwell", "You may redistribute this graphic under the terms of the CC-by-SA 4.0 license.", sep = "\n")) - # write.csv(simd_percents_only, "derivedData/simd_percents_only.csv", row.names=FALSE) ``` @@ -934,7 +935,6 @@ if (file.exists("data/SSSI_SCOTLAND.shp") == FALSE) { unzip("data/SSSI_SCOTLAND_ESRI.zip", exdir = "data") } sssi <- st_read("data/SSSI_SCOTLAND.shp") -sssi_sp <- readOGR("./data", "SSSI_SCOTLAND") # Create simplified version of geometry for visualisations, with a tolerance of 100m sssi_simplified <- st_simplify(sssi, preserveTopology = TRUE, dTolerance = 100) # compare memory usage of objects @@ -1069,7 +1069,54 @@ ecs_scenicareasbeyond500m <- !(st_within(ecs, scenicareas_buf500)) # See: https://nceas.github.io/oss-lessons/parallel-computing-in-r/parallel-computing-in-r.html # Generate dataframe based on SSSI buffers -# TODO: add pubs and grocery store data (using PointX) + +# Load in POI tables so that we can calculate coincidence of other features + +# Sadly, this dataset is paywalled by digimap, though most academic researchers should be able to aquire it through +# an institutional subscription from Digimap, Ordnance Survey Data Download. The POI set is located under +# "Boundary and Location Data" under "Points of Interest," "select visible area" for whole UK then download. +# +# I'm going to take some shortcuts here, preventing this code from being truly reproducible +# because the poi dataset is really huge - I've already downloaded and filtered to specific codes for features +# so code immediately below is untested. Will need to adjust below to the actual code I'm working with eventually. + +## Skip ==> +# poi <- read.csv("data/poi_csv", header = FALSE, sep = "|") +# coordinates(poi) <- c("Feature Easting", "Feature Northing") +# proj4string(poi) <- proj4string(CRS(bng)) +## transform CRS on churches (from BNG) to wgs84 for comparisons below +# poi_wgs <- spTransform(poi, CRS(wgs84)) + +## subsets +#poi_pubs <- poi[!is.na(poi[3] = "01020034")] +#poi_chequecashing <- poi[!is.na(poi[3] = "02090142")] +#poi_pawnbrokers <- poi[!is.na(poi[3] = "02090151")] +#poi_worship <- poi[!is.na(poi[3] = "06340459")] +# ==> to here for now + +# Load in retail data from geolytics dataset +# from here: https://geolytix.co.uk/?retail_points +poi_grocery_wgs <- read.csv("data/retailpoints_version11_dec17.txt", sep = "\t") +# select useful columns +poi_grocery_wgs <- select(poi_grocery_wgs, retailer, store_name, long_wgs, lat_wgs) +# convert to sf +poi_grocery_wgs <- st_as_sf(poi_grocery_wgs, coords = c("long_wgs", "lat_wgs")) +st_crs(poi_grocery_wgs) <- 4326 +poi_grocery <- poi_grocery_wgs %>% st_transform(27700) +# filter out non-Scottish data + poi_grocery_scotland <- st_intersection(poi_grocery, admin_lev1) + +# Load in British pubs from Ordnance survey dataset +poi_pubs <- read.csv("data/poi_pubs.csv", header = FALSE, sep = "|") +# select useful columns +poi_pubs <- select(poi_pubs, V1, V2, V3, V4, V5) +# rename columns to tidier names +colnames(poi_pubs) <- c("refnum", "name", "code", "x", "y") +# convert to sf +poi_pubs <- st_as_sf(poi_pubs, coords = c("x", "y")) +st_crs(poi_pubs) <- 27700 +# filter out non-Scottish data +poi_pubs <- st_intersection(poi_pubs, admin_lev1) # Calculate incidence of ecs within SSSI and within buffers at 50/500m ecs_sssi_row <- c(sum(apply(st_within(ecs, sssi_simplified, sparse=FALSE), 1, any)), sum(apply(st_within(ecs, sssi_buf50, sparse=FALSE), 1, any)), sum(apply(st_within(ecs, sssi_buf500, sparse=FALSE), 1, any))) @@ -1188,6 +1235,7 @@ colnames(forestinv_counts) <- c("Within Woodlands", "...50m", "...500m") # Generate dataframe from rows based on percentages of totals # TODO: fix error generated by ecs_forestinv_row_pct using ecs as sf(). +# Shift to using table summary function ecs_forestinv_row_pct <- ecs_forestinv_row/length(ecs) pow_forestinv_row_pct <- pow_forestinv_row/length(pow_pointX) dtas_forestinv_row_pct <- dtas_forestinv_row/length(dtas) @@ -1369,7 +1417,8 @@ write.csv(dtas, "derivedData/dtas.csv", row.names=FALSE) write.csv(simd, "derivedData/simd.csv", row.names=FALSE) # Output mmd tables using kable -admin_lev1_df <- as_data_frame(admin_lev1[,c(3,5,7,11,13)]) +admin_lev1_df <- st_set_geometry(admin_lev1) +admin_lev1_df <- admin_lev1_df[,c(3,5,7,11,13)] admin_lev1_df %>% kable(format = "html", col.names = colnames(admin_lev1_df)) %>% @@ -1381,8 +1430,7 @@ admin_lev1_df %>% ```{r admin_table_withproportions} -admin_lev1_prop_df <- as_data_frame(admin_lev1[,c(3,5:14,22:25)]) - +admin_lev1_prop_df <- admin_lev1_df[,c(3,5:14,22:25)] admin_lev1_prop_df %>% kable(format = "html", col.names = colnames(admin_lev1_prop_df)) %>% kable_styling(bootstrap_options = c("striped", "hover", "condensed")) %>% @@ -1396,8 +1444,8 @@ admin_lev1_prop_df %>% write.csv(urbanrural, "derivedData/urbanrural.csv", row.names=FALSE) # Output mmd tables using kable - -urbanrural_table <- as_data_frame(urbanrural[,c(1,7:16)]) +urbanrural_table <- st_set_geometry(urbanrural) +urbanrural_table <- urbanrural_table[,c(1,7:16)] urbanrural_table %>% kable(format = "html", col.names = colnames(urbanrural_table)) %>% diff --git a/sbatch_wild.sh b/sbatch_wild.sh deleted file mode 100644 index 58b2511..0000000 --- a/sbatch_wild.sh +++ /dev/null @@ -1,18 +0,0 @@ -#!/bin/bash -#SBATCH --mail-type ALL -#SBATCH --cpus-per-task 1 -#SBATCH --time 180:0 -#SBATCH --ntasks 20 -#SBATCH --qos bbdefault - -module purge; module load bluebear -module load R/3.5.0-iomkl-2018a-X11-20180131 -module load Pandoc/2.3.1 -module load plotly/4.7.1-iomkl-2018a-R-3.5.0 -module load rgeos/0.3-28-iomkl-2018a-R-3.5.0 -module load tmap/2.2-iomkl-2018a-R-3.5.0 - -Rscript -e 'library(rmarkdown); rmarkdown::render("mapping_draft-hpc_optimised_wilderness.Rmd", "html_document")' -# cp mapping_draft.html /rds/projects/2016/kidwellj-01/mapping_environmental_action/ -# cp figures/* /rds/projects/2016/kidwellj-01/mapping_environmental_action/figures -# cp derivedData/* /rds/projects/2016/kidwellj-01/mapping_environmental_action/derivedData diff --git a/wilderness_layers.qgz b/wilderness_layers.qgz deleted file mode 100644 index a9aee007d76683d4598aa3dd43d26041eba10d76..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 11202 zcma)iWl$Z_vM%lt+}$m>1Sb$|;~L!Ef_n(T9X1ZZ-8Sy-Y}__Za0~7ZFXz@d_r5ys z+^V;x>ziIZT|M)or@Gf#p9%mL4i^dv3JEI0%U5Zu!=j^y844;M2?h%9uhrAW-rT~) z(Zbc$$lk=;!o~HYvz6TM_m;dT<49(QB&sK{JxyS3kd(rjSrpGteyRFm$e*e+S9|ZzP3}L?V z+GnGP=Y3%FYVMf6ZWH!&0W-grX3v}AxucCRGrv>X`z)pxSbpCdM{HUezxngp?%h~! zKp(9OeeDh%WwYd)@DE(3y&lX->wVe>T+FOx#zz+YEw4c42utU>Z~RpYa=3|r7A#2% z7%SW@lb5950PuW)5jdycJesEZ`glan{z#>9Q(GX=2+`$W(J4%Av5wW?b6sHcWc3De zh}7Hh%afF~icg?SfsNlG+>H9cBCTIXl;W+5<}Ii^;~7zQp-qD5WHD0-A~9ZUphsi@ZI3g1=TyR6j?^Os z428^lL=MJ%E{ZSNi~k=5zY~hHyv|t9xK5|t`*ZMWGsB5#b-x1LU7tBz;T@72iP$Dq zodq(q`sSi&1Kx&;;_gG!ekQPrIMg!J`fxcmDtSv~F=p$-x%>7vd~X%0<^!RKxeBa> znDY`}@AidG7pJ$ZMOze#xS!7P#$vr2i3l+;A#r4w_Wt-Lne)$3tUJFLouh7lG0ck*u@`yv?4?1r&Wrqqni>@!7RVN=b=P@VV*Gk$E2oltGkMlPyI9y_ z5J`2FdRY{vEH~M1Iin}uLe-#TyO%XflW7~s$}1yAqrX{mu)$v4#X!#!t%ue61JM)HHv!I#L&NHr1PGOE8d_k~8)I?f3l%3M;Pisk0DTO14)d z8VxmnDBsPky4BYj&KZW$hQrKZFlyLd`+C<*|{UY&?~ z(H&iHSg@4S!+h0+_d%vsV||Y7;rJuZFm|9=jmiFr2M3}4xQ*lyzLdlfIzSN1$btys4Lo{SQxqNQS{`e9M{1wU>3icB+m+5b2!JTPDK=Qqyd2)zg1}&+Hmzik~(ivhR+)w-_ zedl;^ekIQ=RBVH2?*L{`;KVPi`r(d_ms3D9pwSp0tXag5*;Z$JyOQ$!>(bOLq&a=l zsS?}4l8vNdz%{5PMsqVn8XNdw1;G7QC2(b>=Mz+=Q9)cu)NW&b#w&9+5HlpFU+;f z^gu^rS6mQy271e*O;&^FP&ZhO^4Umz5tHs7N>POK3{wZgEW$EjJOp2#jZ!- zX&~`Ch0sQvIqgxg#KJxvQCXWg6#$~Kt`?#1WjGRoT6hR7v@`fp^;TM~0Ji9B z_ATdiuf$PH_;I%S^vp>9pg4LK!39-)JLv4Fd!z-GPL4HLOpuQPn#?n&Pa#car%i)@ z^>*hMNR!2cXzxGlvTuLwwS4bD#l>$=k3PC~S(;a=_ALiWYk%}qo@Uru+-hIxR_k7O z|J*kGxk`Q_a2T|zJnS&_b?iEP2gW>q9-}3hT|JTJBYI~E7eeoJ8oUwiF8Gxq?y4ch z7Xl@?5C~UwkBqscc<%I*vr3TxV`xs5!|!R+kHSkyusY1!%_A@Xo*j#`xB@)P0lmWQ zsAzUJgGf3&N!E+yah8>a;(ai%R8oP^%P5;YX+oXkBKQU|zmZ}$&gU7DgCkfp;gE~; zA&WcvoOZSz$3%3G;77-`?bdWnA{$;oz7|Y#_x*8fAJB(6ZEMniSiZHbs!#VTn5w|^ z)XRj4MsHl4D|lGmdyCqrZ5j(04GLjh&woEkDF?Ud8@7*Vv1{P%x4RJFPRYL~niElJ z79WvSQB(sK&E;le$N}o^CSZD z(TVBaDjDXuVh>6?mxj@$6VSfncWOp4ZCa8MURB-GD_$i@`VsW~3P1}PlW1?JZ19t_ zKFXi=8T(A4OX#ydJMD@eM*I;I!unxmb?q~?fJFH?#kcJbko_fCLhN1=M2_eWuoLVp znw$*Op5ARCdTS3y3TKv^3vuD1C1Oz}1Chwgzh1V^CW8j+$vw={sC|)zM+5McJHG*w z&U!>_MPA9Bq;9?wu7670#BADMIZyrwSSPJ}6!0ECVR(O`LKC(jowfVgW ze4`NRMLMhReo163tX6pcmeqQp_NvP9h=1QgA{lL_hWpshCaKgh9p)T`oASE5!ptlb zt*Z?bN<=Z6n9wY*fsE_2ur}5}m{GIpXLNBYdHV3~v@HV1s;X;WA-a^O6TtWEIOxdl zhYugj6lL4KDDYQ+8JlFx;?C^}gwl4!-bTLHSR+|qQ3NFYl(y9M5=N)7;@JRxmcg@HBE#C0@Gv+>1IRRXxde%U!9$puX4NImLf~xc^%f`9Rh; z2Yde)E8jB-Tr9?XJo1R{*}%>pkp z3Ob@PqT=2vJgG}I24l6gb$nYS=A@fbLUC^XH5XS3*WU&%wDwI@miodx%6;o`&Y-D7 zbgXP{R(M;>Mgi!=E)ijuw0dr!IUwTJ@0qnPR}$~qi+tNtgxueeGnt3%mnP?v!7Lj` zQFFajfZwD0TbY4BM{omapCHcQc4g5F)r@;S z`v5xGnT;=93)Ydn7E~I~&Qt-lFaFvpt zy`mb-$6tkZM8W{gfYZDCQ|iP3#GU0djR<%Qo7YGLXR>n=9sjI{^l_@!xPU%_Y&PQT zZ8-sY0rbB91)2KW+1Ub@@L7t_hu^Y}ZNo5}2g01W3;H|HHao7Jr&|LLZ77(Y5v`3@ z7j3h=#1<!RN>3icxS|(}$yiW{(^EnQ;cxX3&71Nbj^)!B|bQd7m z&mAm1cXX4sbFYGYF&8pQyFNHVM8!u$y`F-fb;Em2C)&QFnvx^eh;i;=F?X;-+bR-P zlGkrl+7Rfq-$$m#8C6Zx%>4|g0@lZ$mfO~qu>NR%sN8!pfsC;1CyEZEVkD*s$N z%<2_}eljz}eJ^fs)?-hGF=~wlXc0MVX^pG&fkiJYL8n z!irV!>%oJ!_7_gQ4TRaJC|5aXI+@jW5-1bM;@$A81UIy5ch0qo6+;(ny31-qSqJPN z|I}h^OqZJM>RfR7XOgw8i-#6u$vRI-k)*{9CROKTt=v~1jR{{cRvlwKo>yY2XSAz1 zw7_>yDB8&6?h<`Fi?I|HolLRaW~V^tq*I=1Ejgk$(I1V-Zx)<%6ZC1z@GI64wYOs5 zXr!$mX$`s0p7L7_P51f1_fjTfO+|hGfz3QgSQ_ zakkyvx6AQ$Ub`prOV~#QIT38n!;E!Hoa)xz2>jYmU)5KqVguO2b6!Y6VMVmjW9EswXwopm zCwD9p5N%phOYk#pE218gI88xegW0aa8-br&)zEHO-`+jNFR^bHXmoIPJ2n1JxNrL{ z1@nwIhzHjwgK4_zu=P=y?a_hTh3_-p_DE1(_2{@0`>KWGFZfc!Z6-!uJ_t({_%w8y=SEe1bdH zy)KMQD7*wJ_mC#HUah)^klnsNb}K2xImNyo8wn!5pXc^tfB)@a3BKr|_56KudfNH; z`a1Z7wlsPt_CP`KeiHUfAnUkfmwc#c%f1hxb@#WNAT$+GOwERO$M7zBarC!$bh*dj z@Ang#Q=FlQ;yN{fL* zH+GZB>>MK~BEW?r`bDgW?$q;6RM-FkN8z6B{bj%V{n=WTOkTidMWLRFte^Xr+Yxy+ zv;~yU@TA^OOZqLzipr)i7w&>_I@1g-}X3ql|>Qq~lF6WmP87*0smsBpPl+%IMvr6|$S2oM&Kr5e3 zwE9S-a(*^f)AMh+NmyT&;*#Aj&~zB2)-0&~Dch5tm#499#B%TWc#EA`R%m3xn=A}C zDGvRfzNb;g+L$FDOMF~bT-BH(e}#vy0akaCOx}fPHCXe!4)X{_7|~tBdPu`54lo`l ziq!d{Hm0>JmRv$iQwl%y`@68JaI%Y1vrB$^yS{HN+{rZdk-`ZE$MJi&aRjK%fwHc+ zecEW$^U!}?2Yny^B-TbxQ^sN;b8Vecc8H*6hw-LruE^U~3|=sbF9b&x*_Q5aQjmc& zW_eYb3~2!PHOIV&?>*|~X$FduBMD_=ZEZunV(hobVWZrEx_-JlMVRL{x(v>sIDti( zU&l0)w^#>H;jaDzt)lJ+PvQ*WbABEtJQNL+I->o?Ji0ryIj*1YPF(%;pL9tCDno?A zq_yv{)b12rDQrfECQp|061{L8a$et~COG_jdj!xP~p=A1k zbpP`ceePA2%k>H##^I@)=$z+agintQb_YEKtxM@)*6~?W>DEDTujA6hy9ciFE0jVx zb3pB!_^f-ED71-0^j%&n-RaM8C3tLw=DM?O)Eagx04m)Zwbts82q&hwL-We8>4&iX1M5nQ%yjEn+B_F9UJQ zO?fC(1xZ4ZN>W-}vb6vP+<3YY-5h2*U&Q$NZL-)^HkOii8PH9xZqW_px(dN8N2z>I z*s&RwMnjy77ea`fQyE%4`aEa z1Mu?t{||V=(1CLOX#vk8NpSz+F}(cCpz8JkK@I}{-S7^kl@Ofz>*>kv{=?j|v9_qU zhqt?npAAJo*P}yf&x3Y}Rr_pD&_-_G{-xCAvny7j&ZhFGHMTLn_eejJ6)r!>;o>etsVU{NdqW7WlGW;Bk;UI^m}5^4 z>uGOhlNd$aXeom69&W3$NBxuB95QpD|D+XVzSEl$xr>;1+O@~0F8^S1qPeGc5t}a; zDq(V}*y5GwkhCFj3Q|RYvl5XZ zD762=v6VRGVirhT`9W!ERt0w{-5f#OsQhzlLd`rYyb>qD{Yt)034gz?F^QLZQ1>ZA zCxFd=@1-lH&63_IM@CKL!0?MvMii*nC@rQspM4^YKYmedE-SeKAk+Xfvq{x*VLh6c ze=1tn%_^y@G)r~X@I3vSG?NB@lcrNBRYcvbWh7_0^(7s4d6rl1TsJbqzr1%Ey-?4c z!ZzM~@KCor!ynK)&4$k8{)s1)ssKc9krT+lw%mYbCx?OlXc?Ryup}|3(66>Ql1GI0 z@uU|Fnb^A=R@X~5SfBHs*QW3b9Nsd=Rf zhr;L62oz(MYYWqk&3t2$VzfiPb${x5>`?cHBjC5Xb85e-4f}guVK5txnL9u+<>*m5 zbUp<-8LLRl5&KeL<08EuOT24Y_?8O}!&VEjkI#@`%a>4-E(W4=adBmN-7;Z2B+}H0 zItF!F`Tv&wTn@`7rIKacn{y( zEM_|)s`R@iVBdUea=3oCo*Wp@uxAXuL%(Ot2)VR~_MV#SoV^nv`d(B%`H6SJ9!`$O3QvTJo9NW|Avef28n|Q4N*`>e2va zX)7wyiCXh)ftH}l)}QPIsU^4;S^ZV$h}{#JNSUw2uj*rq~9(}>UXp; z4Bt3qH4pm`C?S&943@z!R|}d~Z4V9hl1JApaKk7BmFeugmtzUG&GYhfx#{q01-ZTqL^wcF z?2&DjAo>%0m(vZJ%1!o&?Z7G~$>n%7qhwNzFfi55>MCwpKb%Jqh(PQ`$id{{vU>CM z{q*@$3|#kY6oZE|fPG}a1{yi9c~()YiIp39xig()r_k+qMy})8m66=BgXf9Cw_Y*u zyY$wsBb?c7waV%8!^PNhU>PW|T~$L7G`LZvYtdS{hH`{E32Fj~m zM>9y2;)VucZ4&+HjBAGvlhRGZh0T&r;CA0G(MtD#ms&l@0#6N(h)7%U!2keQzq~qi03a z6^7=$mjLatB65?FH@6c45Lwa`R-8VHY$442yMv4Qp+jQ-XYQ$`&5r4LUS5HZJoo4< z-qF9ptS97M5*bpp=k(x-s|&NM$^YmoUE7W-LaN~37x~7U9n!dO?OlDC)fe%xU2Ywy zf~<9fJX-O-wiceT!|Tv1eFHF>?n^XlEG{1{lPXoJ>I}>>s%s z1vZB*vd1M=qckd$D&e4sUK1+{bn(%2J4k1fm=x?ds{?HN)YBw!KCWz`6YYlTU|F1E zw_=fu#ZWfQFxLheZTc05SF=4j>0KS%LC&N$<~MbbnOJL*32o^d=tdmW4M(L}u#Kd) zqh$TO*&^5>+?^Jznu#a6+oRxbuX0^ciRZd2jmPaAM4+OVHuB|DAM@S}MGZs8WvPtm zNYN_u@^gs%72T_@zuz~5DJ?brQcuNmt=aDesTUpwea5EOyEcu$x zZTYtHiyRh!dREOkT)BE41CS7gj=Fih9btJ9Orr_&h5Bbp_xxm2X8O9e40@vCKsQC8 z#y0q;gE#h2xvF)L(x*Kh^Gw~I!jM;EdI#Kkf$4N*8+j4PC-p`SFbdt=; zX)TK{jcgm94(j{{r9ekxupYy+t%M+$;`vdnmf1+1bchwh3_pxB%8xUFo&IW2sw-`T z0O!%Pn$1Wdg$>E)&O(Ra(UjUqVK*w`-(0ATPD&;SaBLLWlU#8W{_aq?bK4OxeP-8K zl{s0VonOsk7V1*g68l>!JUO5vPjz@!S_RM3Y@>gfr`n=Qo?t-uzW6SIr*bmx8|5s2 zJ0vy;3%>DC;7S__{9@P zQeV;NM!*WCS`__am~LBD)2TP7;w^40Z%nZuBzjy!H=I$Sq%pq9Nm8K#-Pc=dIjxhU zBUN7Zg$U<*8Ym<78 zxL`MMUX;qM-X>eqPPlP^7nPIc@IX7j7JQ&>{RjDA!SoL@ipfillMVn~Cb zMC({#Rby_VByQnYiSAMDRHsk1N-P--Lad>vTP6u99BU{A7RbkWe6+Yb7Bd_=IXsHp zjO2^wj{2@JSy`Od7WaV}Sg8wHt|_Ntiqy|~jx^U%qn(FnF|-5c#}l|OPq7pw@E}9zt!YhnsmflI$TIJgZWYQ=!e8p^^stp_2NkXBn1l3*l!&_j9R}qE$?~yE&0B z>InhNSIS@B%0|TS9!uOGOQw;jk8iBT1<4Q!W{5wP90Z>%6P0}QuR^EAzoOvOxJzI< zrwP{DRp1-@_$o6(6w|M|$qyho9epak7D+9siCmLSAYpQ$`QlhO5`){CP~}3CFcP!o z=eQf|pBHQY`Ao(Ak`=|Xa&%t%|Ac=5+y&w`6I&Pr^+KEe|*{1sZTSO zHoZCg4i)pTsX4-qg!kQ$YFJLTRW*|kT_8NVB&Xy10c{Qye5kRi*NAlNPOzHXgI-mzw z4+1q>zHxl;A`9aUeVHC?@T~CDlsRLj5m&u_DWJUYYlbd5>>QnD;vNVzMpVmc$4eon zqKORFB47#yTchh}nyXhpM^YJlQNe;gVnvYzGBHImS0lvs)^M2|q6^9BftyTOijRFl zOO|DM{q!^y#`x$K8{LP+iBsGmkYDayi7F#vjJux6Ip+XmUUudhTbK&^{~NaTkpD{c zsU9-%N?47zk>-5@ZIm)y_`d{6Aom9>%0J;jmRzY~Lo09ap9q39E-b5z6 zNPMM3^N*Y$*`*XgCQ?xb%(}xbw;iN4MTSbQ4gxIzvMCB*T+tx)%%e80Cx#&{y}m*f zk3Ic~Jv3uCwK1%yx@!MA@1Or5ROir2{8I3hxuFp@U>zmFiUF2Hn(3KKM3%b=Pw(e~ z8(Fl%#9cz9JtS3QRg-sWD3p$BBc5!m%Ev$~aHX_|oeI4@1Jx-cOoiy*`=;Iv1bba7 z$H;=(9f>7WSY$jTUTz0_9cC`~{-I)SI7JPORq8gYKgAU`|ELjWeY5CU`gsKSsn3fa zXmIhpPKG`}(Fvd@ti?mosZpxlv2xSJ)(2#nPJ(Do;OT^fv+1@7p9WbX6!M7@(nMHN zeh2!W)T}28>`@BU$hdTC`=n^+w{5=NURr#_Ek^z>zJO@T1-eq>CHy$19F-w8xx;XQ z4mWQAA%3}z1kTQhKK3Q4nJ`Jz4u@uqv|F`*5hmccj6->Pj{nhv8mO79PnlsQMjrE# zN!kzDQ%U_X!~ASmli+AC6-{6X*QSmol(m#GscZmd!Y`&eD;#QSTNgRmt+mbEcIF+* z5i#0!e5Wm({Yjyw+clGrBXT#c#UNk6nVuYE!YDEsX~f7=A^q`|HRxdi@N)Itxct@d z%}>eP2i#=b+6lx8o;=SM6__i>*eAEz_JV$ehJEz!>xV|1?wkya(y^3~gU(GLK;k;W zgLyu;g0B9`kyrAZd~PhS{8Y0G#}i|_D2#3;b(Eq}g=DLwDaFh6TZT}X`B+s_nv7ALYR@p~Cvklmar!qhFQ|KW zZh~!n(ETyFF`7>L@a2RR@<@+%7%#?9%XgI!8@60kyTe&B{IZiwZUtaGc386ZYBriJ z@)oDF7#>4!VQyTZ_m_9=`_Qu#TJl_xSvGTwv!2)9Z9k22D8XpO+S;T_a_hlzlk{-q zAdlL^f$hDJOk7DV2BCwK52!$>pY##GEnzq=!P&PEm74HQe`8DAf*GlSF&$YWb2iFq zSY;GEfju&{$ni=Sd&}s36FmBU>Wv#a)s-01mll+MoB>J`9*s|4A8Bf>2U@I2)MfcCL?( zstB*i@_F6qzu+l=w&TeFgiRD~5esH-wrC#SSv>T*AER0Qv6}pz9yx>`JyS!;2E(Ts zQrX~L0SOS3-lq#=hCqpT97W+zVc#H3xxMOL&MC&;aS_27nc-RW=s{^?CYN($b~P%8 zL3Rw#TqRx1-0qbc8GYj)GK?a01@pks$GI$*61Z4kSHLM3-8r)t7?z)BL7vr0%#8-4Y6DAlw+F>i3u$!ogmuixaTn)O7x(7zTry1O?-~t^f3QQQdz6 zK~-&`REi+UUGj{jkHi$^1&B*)b6L1kfDLe{61Do{e#tboWQH0jVD!)E&xNeU*$As2 zX)Qz}9ITCI=#Q<-Lyqm|2PqBCMt7qqif18I2w&$QNyO~zcrb)rDoQgeUImJe`MSyt zu@t(eAuA8Lw<6Q@FUraHI{o*o*TRk852BUO1b3m>xU#7X&dUW~pe?lU=8_WZ|M+H| z-&_V-mdSFpM+PfM(pZ`IZdy?SmXG7P>ttav%$XzSA(|6#mW4zkJc8f%Ija2>wk~`jq+-)0)U3Wh5g@ovj2+E s6!;JAKeA^3C&Yh#@Bf6j{!769Pa*z(1w6t(++hA*p?~l6EdQAOALkCD)c^nh