From 644d591fc89b54f1177a785e8f17aaf8d61a1c14 Mon Sep 17 00:00:00 2001 From: Ben Varick Date: Wed, 3 Apr 2024 13:21:15 -0500 Subject: [PATCH] added crash_summary_charts.R script --- figures/crash_summaries/counties_year.pdf | Bin 0 -> 7490 bytes .../dynamic_crash_maps/dynamic_crash_map.html | 6 +- scripts/crash_summary_charts.R | 130 ++++++++++++++++++ scripts/dynamic_crash_map.R | 31 +++-- 4 files changed, 150 insertions(+), 17 deletions(-) create mode 100644 figures/crash_summaries/counties_year.pdf create mode 100644 scripts/crash_summary_charts.R diff --git a/figures/crash_summaries/counties_year.pdf b/figures/crash_summaries/counties_year.pdf new file mode 100644 index 0000000000000000000000000000000000000000..ef1418f56a3a6cc51224b7ccc9ba7091a81f7686 GIT binary patch literal 7490 zcmZ{pcUV)~^7knMiXzfPY66IW0wJMG?+DV7-eL$v0)Z5%p-Trr5k#5@Qbj;OdJ`#1 z?^2{ll_I@_`U~FUx#!;dzWbkK&6>4l);!N%dp^&XTUS|C04yj%$sO<`U^bv9V88-J z2?h!Qkq;kHN=s9Mlu&Re#u4eF48_2KyvlbWLJ$!l5g}o)FyWsVA0a$57{Uoast7`Hl;JR>E&M!h zG=@+!XG)KhwPbT7TN{=C)b5?Vyx0IQ;r)(+=O5)GM4aD!;tiBeX>Tjhx725=QaHtV zrp{=55NHeWC-L35Ckt{9Y38<9|1}1(2KnMW@jZQvq3fqx2uT@M{E%eDx2fiz^Ks@A z`rbrV5jkfT>xsr{XZ}+2JwG;(KNn-BBj>kpbQ|-a^eHqRpXk@zP}otVHg!hiRV%&K zF|x6^$fYxYL+|Xcq2*X8TPw3a&egV-4YC?1CakWQov<-{A5$_}1?3-fA6eHS$y~bJ z-Zy?XBYygQ!fYGmW@v7Tz9uc_C77sm(P?|TdV@ZUyRcB3c!n?5($xDE$tT|8a-qnh z)seC6%HqDW^g^ z25Qkakxo?v@@ZIQNHSGTltdq8RVO?eid`9ZB}-PXjF3C5waZl1W~=kajf0w))#cLEmK`MRTPqd4*6U2_C;c~MJpD9JF#S}1dLrZqr=U;3Ra{d`8qL5wVYUAN$vLcq&A zw@XVx_CJK@i8iX#ZBnf}t}s0ksklUcBg>52hNelH_M-481d5-d-5v?KXO2t18yv+v zE);%4wN-9lY^6}>p4fd>Uc0HD@_=m7}=VmI&Wvzj5op?;}4IwItnj_{F+SkDdy_-wz zdrm8yG~ke^U3!YoUr)PDn&2AmfamlIE^lt4zUA+uD54cl`T}AX&@dBh zd!G}3$^WN;x(Tzj(dOA%dPG`2s_(YQd#PB=hJW6BSqPU3LR@r$m#u zL)*(@BN;N|pm_ZReq5HVyq9Y#@uBV;Ka z!!swBTXlYyy&peRTHg~?U3)L5G*q5aYj(I~(MJt^lK&wi^>+W4maUV`(v<=pT_;dm zgOw~EDv3W_(&5>U-%M5=TOXHYxgNnm@8~cnKXgy`DUo0cDocQ(>Y`eLj^2&0E^{eQ zyt&>hNV6|n&ZPzC^*s>%MB8Ira>d>yrt;@WnW#8__02cZo3RZ~Ve^W+ve}9+@@pgk zt}&{B>4fNd*M|qT!=(IDw~81mFTkj$$D-@}y@mFae?3GlTS<(9=kX@lpG@ZiI6igm zYJ{4{XJ3J=I<#v^7W>lbkx?jiD^zZIDKd0Rr$t#W<$ks@WW0(<&y~{O_^jQQ*=)Qf zAf5KKoiwp}{zliq#<3^G3w(g4*dcKnVu#}K&^zkmLX}eAY1(L9uU~Qr;Lu%o@4=a0 zEn7zEl)(tm6()opdN2N@4QZK z&Wx+P+)Sq`=v4VryN z%uTflu}gopS=@f0PufROX^z7^A5ioGyZM+ZP=Mj$slGu{c zpU#g(OAE1HOEm4hAB>Y@`#H7Cw)vu zM2m>k#jG#ye-|VgXafMtmk_xr<~3QT^SchSu7irM7P)gtR*nyy9>{RmVI~ci)>_H% z{L%>M9`su|l>v?FJ21Uz+JJBrA3G@5$KEd8kvNDtN!)$<^9c6>Nl)TUQ-__qSCXhf zbK3mi>%A2HynxA1CcASLLNTF|3gX*PBS5sQsE+3#V_zWrgg!ao?nkV>>Md1Um52eOzlbS=xv< z=UBM$>2Qcfvr~0KIu73Rx|U^)58Ux!x46P+^w39{y(8!&d~(3n$7AwsGClwXK7`_v zqEAiIj7sZixksPxl6IuE?J1nCwfZhi4c9N9J^jcdFE-z&O*R|67mV!g_M7*D?k`N~ z*wpx+K#`Q8<3+Dzn0OR?Ct)1c8Gf_vHVRM@xy7}|I~*5!jGIvNMJ2k`zJC(iNnw%Nm*eCpcu1m#93If} z?Z+RsKcz$$wOZAMuZPP7)Gg4xC_1`}n{rib*21eR6vwps= zNd0(oW%5bqB6l$?VN>Qsr43**xy0HijC+swOfkXax$mnna}Q}i9CGYYQm1zw;i8Eo zTr?RfnE;`5wl{LF>UOwBi-YzRPHu3^8neVNa1mxN zkQ&n$SP5X21w7X&rX!Zc~CcdR5`)o7D*FY zIaf$0>iQ}qoV_=^{0+b>0@^H*2x`LJLbynT_xP;}C_c(T66cS7_+fi8b0u!qSFw-o z(>i5RvnDJUgg6tMGDi+psV1~O9Nblhf%QkJ_-EB7qI=oab-2=tGE#a+0 z|MmLtN4GQ)<~9Q?A=^ABMk+g`+LKG^81eXB?(ZkIC->atN~j;xJtJHvXvpukXF1@| zmBZVSpWUxQN~SbAc--33=h++L);z*jpw{0li;V6|)3tJ7*v(Z@3ZiY2F%+@T&v0h- z$}G=DlE>&%BdzZT9^{=1JLHI4n?>30D5HVvdG!Byg}(D|s(gvbdOj zxGX^og)!kP%yoK2Of!C8HRhv6DbQv%%8TezgNlF-Z+i%I-ub00*>&;4QT*>&@ zLhS^C&1{I*s_ibUpuC%bctL37a(g(ZIVtXgmU%`Svu3`WoE@3GQs%k=@kU@WiSSWH zT|f;$3>PRwJ^0ctIYxrIxs`P(H5ZEz3!}*m1blF_&?d@vMNL zk!=5d>_T-;tNWc-JzkHZ#Mc@~W-pDZ`gP^iD=_J(zpKa>+Ph;uoZ6_@3SJiVNwf@m zE9O3UfIQ6CnZ#0iuhbxMvYghXzbJy>O4FNos+APL{!_$GGu2I|lcw~n^q+%gw#L)D zvk8BHH!8{rU)tHVD{+V2FM(T;L5%i^iAMl;Im@3J$od(-jr2mJzM1h*{E<~*alw@q zhV0qFVu={-@}{^cKc}UuC+x^BfQM5AeYl4cw!+zFuTg%#d#dJ)a;97Scn&g~jPK&v z*^#+h6YmRYzz6zT0jxWbvab`V6#zSnnb$d{WGLS82#QTDjX&L*-&$b79ck5DO!>S) z%6&;(Ne1n#G2&6U?ST3dp!IuXeNV;IY{zcbB7Qg+@U!8dDJV!S6EqvRz1w7E^3!s8 zVPUT~nbcNE)*1PWx=!+=hBaSOa&u4+5BBu+_sLa&7y4JY<@f|EyBGRr^jPuU!U0F} z2n8Zz63Ji_Q7z`$5f0+=G*tpPj5aU4yb)HSeQ_guC2HwxP`K z*itEL?VvJ-Eec*u6I{*oQ9&JA<896CJ-VPIMlm{%S_@l|Lw~g6j8+W)7rm?3Wck!C}q1~aDtscgX zOZIzHrzdz+{gK=K)5(JtzZzmuYO#2S-jx%JrPD98(EH`x9o)XePjB^kS*Ka9>ma+| z3nmsAcXTk|*+TGgdV9@qQiAyn0<9{*elr`Yj-fOa- zPZ;72XEts_(s z=>aqsBFu?PLV#jmh~+;>_MAljkS76{VFqvv&>TcSc0vTW2gZ^Tq)MPzptLkd6#@jG zZ}OXiK`I^?H3Q5!u}e#zZ$}7;K%xv>p|IZ;5L#bN5h(h{{omw#ZaR-=4!WnLqzFaB zZGmFvUY37Y{CxE1j0)%Q>vA3hqzr%T2!o$*3eqOXiuGZo$j-15%IlF4*DeORCVfcg15n!u4y#J5^;(t8@ej%#bsZF@J{<@D#ev> zMm!6O3#vT&$f_hs@8M8D|rgA@9UIa*Yn+e9XDJ}6nh**_UXRKCI2fK zCiF9lEn(4-_BIs6S2l=95k%9jrKO734^&B#NPNCUA(H4080)I90`W(;2jrrfh|bdE zt4uk;M7l4d7Pk&YEtnX;NXQtca98MjdHGR9t6NB;o6|H9V)Pun#c%GjmjtJl+=c{8 zGj*W5f5129CU|C}C28xU;%t~07cU-ELF#Y8f-v_&pCIlDI2DS(R z*fU;8punsr(X{#(>y3VpNoNmT2X#w#r6Uxoo>B7j86hq_>(qA1=H=kyWRwaaw9nkb9n2@~XV3zk~C*r+L z1H9yId#86Dcez!K6yFuum-I6As3-6gXIy#zwEg-`zZjUBGNh7LuY2{?u^EX7kVr+I zbXlHdiIlGK!pBCMN+RlMQW6E;unTTB;ORj6SYl-$BfWq5>l^Wbg*H67?X5f8U$IXjgG{{ zm|(`HFB=dEfM8%n;|o_(pKHH@b(#~WX-p{P0&-~luhHCn6DdS<{gbLeG#mx7 zv{@vM@?){=S%9ht^|+jNb9!*QXxrP3vu9f?B*@}LY@6gcLF-+k6w45cKcOr|nA0ogXDZQ_G(|VIh#o&UteqdZbN?2P+Ylq$W{9xM5 zmB-cCm9CPhmUi;mdK7ib_D07h{6N)*+V`SdOgaZIUDT!b!Rm^)p#mAtdQ~MfD>Xe7 zZqn|CX@7Qo<6jwKa!*=WT8EXpo?B1d?xq-Dz)cQkjy0ZC&hMOcz~^67ZK+YxruQcm zU2dgwOz|}6y}oU7`vUI+l-<-f^JuNUsD~GR z)-!u+n29m`8PZ<{c`&LfX(*}WS&E9>BAZvE=Whx7p){ix!F$AO&kHjQ&ewq(R1IPW zRhLbdC8(5yv&nk(9p9EA9fHdE7$PE4HOe&7G@`${G4Qcvu=J5G{oH{ZLJzTrQV#?k z=sc)jzBME^bSJ-PXd%}#cXLQ%aBN6*a57JL5H~n5Sdr(RH#H#d_|h@hY3|W#pG3hX z=Bjf{2kWntXOf|d6T`uxT_Ri}#ty#*SbI~3`-is%zV_$$-_6`LOn$$Xk(ytpbQHXo zb@f4;mfGu5rNyjnSO{hF^SJI@+Z}K=MihI3<*weX4yqRRYW33ddg8U_wZC)ekpJZ7 z3EN@ap#k|_at?Ab@_@ONjg(JNx5;+o(BFbgEKi>;Gvp{6TTp!d4=b;U^XNS6Dw1J9G1 z5%3tZet{r?%K|n0>84C3E#}Up%N2(X&d5}k?+&=OkcV^*1`Y$T?SjH2fsv1`dt(+` z7E@IPi;sJ7F7_^M72`1XfwV_}hYd~jBd0+@nL+ItUC{o%MsQ!r@Fsa=D<#a?F4bL8 z%2TrN@!f@;;n613RiAHoM#gN-SDKL>nVzlN2iwex(N)V;D2o-#w=P#Ilq>37^xZAp zW4wcoNU3quX==~gf^+ih^A{C*?5&2aI0?8rg0SDy*yBTf9`C97?A;ssgKpkn4LxPt zBR*O&m168U6`&#%lc_Dt-Ih1iUE!hgX(Prj}3MF37A>ou9ekOFmJSD;Fh42P?_ zP5WMK%{JBIeb@onm`!iUDKRuKbXbf2+pu-VMnV z6v$`mUqOpHH>+;F;E~aItTmY?o%AD3KW@vg@w(yVw4B47;kWhLm$k;z_xNOP)ya@1 z#=UeiRHw<66<&I8%yk&y?AfGBhJ4>B7qrO*fN_Kr$%p(~hR9eKvDiOKkVPm}CBq|+$ zW#?-FG(ET(H48QGXiOm&TnK!9=f z4|Qk?nA(?^M5)BpT3x8IUPWZF<_AqPCW}!^xNdp*G7hV;s~&Bv`ezvP3p-DZVQ@2ABvWjOl- zOJ|Dm%q+TPi%H3$7yaC8F)Qq{%pPSGZEtgCy*t+T0Ru8dXp4D|(MzAg-_WINFQli) zzQ8X3{5dxOtJCyaPYIVn7lLY$8;IV~?M0uS)SvsSa-NgB$0e7$8Dp$$Emy1OT=US4 zzOQE~jsSPb*~q0rSwpKqQkK7%(!_m!L>#>sFs>^~tGVvF*57^{H?yKQf}B*ZOH?G}y)BxiNtaeQz1c~VY`+?_Lvqt3Ze_Xdw+h!BI6=&8ld@jy|s zb#BAAv;JT3b4!(%}aMW;i~nL^*eGeeV6bzb~bCX9p+fK<7e_ zIZ#9tg?2;(L4qzwj4j*_2txf% z5=Q9!5eVo%><(BLdngL)jDTV>KoHU%=>m5mM1%g`8GYWlTmC~6``zpRX~h2}iv1=5 zkOB4~=C?+6ZX_UvHq@Dr5`hkW*X3V6>TtwkIK~kMr3C${3-G+6|J2YtP;fg+As`q` zDfIUP6cZH{76sY?|F%I0O5ynfbotu`79()bf7*m4!341QZ<{1RR{Kwz2>AS1|7{cb zpZ$a&V1g$2pMF9R0zLhgO-w@ge`1IWN&X`j3<~OqfTIXPClF-d=mjU77YMrNVuvKq k*?B!c`bZ=OcrHNxt_>OkMPYtdMjR{&q2%UP(NU%RKZV`BH~;_u literal 0 HcmV?d00001 diff --git a/figures/dynamic_crash_maps/dynamic_crash_map.html b/figures/dynamic_crash_maps/dynamic_crash_map.html index d186c18..471bc03 100644 --- a/figures/dynamic_crash_maps/dynamic_crash_map.html +++ b/figures/dynamic_crash_maps/dynamic_crash_map.html @@ -5254,9 +5254,9 @@ function asArray(value) {
-
+
- - + + diff --git a/scripts/crash_summary_charts.R b/scripts/crash_summary_charts.R new file mode 100644 index 0000000..506d0be --- /dev/null +++ b/scripts/crash_summary_charts.R @@ -0,0 +1,130 @@ +library(tidyverse) +library(RColorBrewer) +library(tidycensus) +library(ggrepel) + +## Load TOPS data ---- +## To load TOPS data for the whole state for crashes involving bikes and pedestrians): +## Step 1 - download csv from the TOPS Data Retrieval Tool with the query: SELECT * FROM DTCRPRD.SUMMARY_COMBINED C WHERE C.CRSHDATE BETWEEN TO_DATE('2023-JAN','YYYY-MM') AND LAST_DAY(TO_DATE('2023-DEC','YYYY-MM')) AND (C.BIKEFLAG = 'Y' OR C.PEDFLAG = 'Y') ORDER BY C.DOCTNMBR +## Step 2 - include RACE1 and RACE2 for download in preferences +## Step 3 - save the csv in the "data" directory as crash-data-download_2023.csv +TOPS_data <- as.list(NULL) +for (file in list.files(path = "data/TOPS/", pattern = "crash-data-download")) { + message(paste("importing data from file: ", file)) + year <- substr(file, 21, 24) + csv_run <- read_csv(file = paste0("data/TOPS/",file), col_types = cols(.default = "c")) + TOPS_data[[file]] <- csv_run +} +rm(csv_run, file, year) +TOPS_data <- bind_rows(TOPS_data) + +## clean up data ---- +TOPS_data <- TOPS_data %>% + mutate(date = mdy(CRSHDATE), + age1 = as.double(AGE1), + age2 = as.double(AGE2), + latitude = as.double(LATDECDG), + longitude = as.double(LONDECDG)) %>% + mutate(month = month(date, label = TRUE), + year = as.factor(year(date))) + +# Injury Severy Index and Color ----- +injury_severity <- data.frame(InjSevName = c("No apparent injury", "Possible Injury", "Suspected Minor Injury","Suspected Serious Injury","Fatality"), + code = c("O", "C", "B", "A", "K"), + color = c("#fafa6e", "#edc346", "#d88d2d", "#bd5721", "#9b1c1c")) + +TOPS_data <- left_join(TOPS_data, injury_severity %>% select(InjSevName, code), join_by(INJSVR1 == code)) %>% + mutate(InjSevName = factor(InjSevName, levels = injury_severity$InjSevName)) %>% + rename(InjSevName1 = InjSevName) +TOPS_data <- left_join(TOPS_data, injury_severity %>% select(InjSevName, code), join_by(INJSVR2 == code)) %>% + mutate(InjSevName = factor(InjSevName, levels = injury_severity$InjSevName)) %>% + rename(InjSevName2 = InjSevName) + +TOPS_data <- TOPS_data %>% mutate(ped_inj = ifelse(ROLE1 %in% c("BIKE", "PED"), + INJSVR1, + ifelse(ROLE2 %in% c("BIKE", "PED"), + INJSVR2, + NA))) + +TOPS_data <- left_join(TOPS_data, injury_severity %>% select(InjSevName, code), join_by(ped_inj == code)) %>% + mutate(InjSevName = factor(InjSevName, levels = injury_severity$InjSevName)) %>% + rename(ped_inj_name = InjSevName) + +# Race names +race <- data.frame(race_name = c("Asian", "Black", "Indian","Hispanic","White"), + code = c("A", "B", "I", "H", "W")) + +TOPS_data <- left_join(TOPS_data, race %>% select(race_name, code), join_by(RACE1 == code)) %>% rename(race_name1 = race_name) +TOPS_data <- left_join(TOPS_data, race %>% select(race_name, code), join_by(RACE2 == code)) %>% rename(race_name2 = race_name) + +## make mutate TOPS_data +TOPS_data <- TOPS_data %>% + mutate(Year = year, + PedestrianInjurySeverity = ped_inj_name, + CrashDate = CRSHDATE, + CrashTime = CRSHTIME, + County = CNTYNAME, + Street = ONSTR, + CrossStreet = ATSTR) %>% + mutate(PedestrianAge = ifelse(ROLE1 %in% c("BIKE", "PED"), age1, age2)) + +# add population census data ---- +census_api_key(key = substr(read_file(file = "api_keys/census_api_key"), 1, 40)) +county_populations <- get_estimates(geography = "county", year = 2022, product = "population", state = "Wisconsin") %>% + filter(variable == "POPESTIMATE") %>% + mutate(County = str_to_upper(str_replace(NAME, " County, Wisconsin", ""))) + + +## generate county charts ---- +county_focus <- unique(TOPS_data %>% + group_by(CNTYNAME) %>% + summarise(TotalCrashes = n()) %>% + slice_max(TotalCrashes, n = 8) %>% + pull(CNTYNAME)) + + +TOPS_data %>% + group_by(CNTYNAME, Year) %>% + summarise(TotalCrashes = n()) %>% + mutate(County = CNTYNAME) %>% + left_join(county_populations, join_by("County")) %>% + mutate(CrashesPerPopulation = TotalCrashes/value*100000) %>% + filter(County %in% county_focus) %>% + ggplot() + + geom_line(aes(x = Year, + y = CrashesPerPopulation, + color = str_to_title(CNTYNAME), + group = CNTYNAME), + size = 1) + + geom_label_repel(data = TOPS_data %>% + group_by(CNTYNAME, Year) %>% + summarise(TotalCrashes = n()) %>% + mutate(County = CNTYNAME) %>% + left_join(county_populations, join_by("County")) %>% + mutate(CrashesPerPopulation = TotalCrashes/value*100000) %>% + filter(County %in% county_focus, + Year == 2023), + aes(x = Year, + y = CrashesPerPopulation, + label = str_to_title(County), + fill = County), + size=3, + min.segment.length=0, + segment.size = 0.25, + nudge_x=0.5, + direction="y") + + scale_color_brewer(type = "qual", guide = NULL) + + scale_fill_brewer(type = "qual", guide = NULL) + + scale_x_discrete(expand = expansion(add = c(0.5,0.75))) + + labs(title = "Drivers crashing into pedestrians & bicyclists per 100,000 residents", + subtitle = "2017-2023", + x = "Year", + y = "Total crashes per year per 100,000 residents", + color = "County", + caption = "data from UW TOPS lab\nretrieved 3/2024 per direction of the WisDOT Bureau of Transportation Safety") + + theme(plot.caption = element_text(color = "grey")) + +ggsave(file = paste0("figures/crash_summaries/counties_year.pdf"), + height = 8.5, + width = 11, + units = "in") diff --git a/scripts/dynamic_crash_map.R b/scripts/dynamic_crash_map.R index 95a6f05..62bced8 100644 --- a/scripts/dynamic_crash_map.R +++ b/scripts/dynamic_crash_map.R @@ -1,6 +1,6 @@ library(tidyverse) library(sf) -library(tmap) +#library(tmap) library(leaflet) library(RColorBrewer) library(tidycensus) @@ -81,17 +81,17 @@ focus_columns <- c("PedestrianInjurySeverity", "CrashDate", "CrashTime", "County focus_county <- "DANE" ## generate map with tmap ---- -tmap_mode("view") - -Pedestrian_Crash_Data <- TOPS_geom %>% -# filter(CNTYNAME == focus_county) %>% - select(all_of(focus_columns)) - -tm_basemap("Stadia.AlidadeSmooth") + - tm_shape(Pedestrian_Crash_Data) + - tm_dots("PedestrianInjurySeverity", palette = injury_severity$color, popup.vars = focus_columns) - -tmap_save(file = "figures/dynamic_crash_maps/dynamic_crash_map.html") +# tmap_mode("view") +# +# Pedestrian_Crash_Data <- TOPS_geom %>% +# # filter(CNTYNAME == focus_county) %>% +# select(all_of(focus_columns)) +# +# tm_basemap("Stadia.AlidadeSmooth") + +# tm_shape(Pedestrian_Crash_Data) + +# tm_dots("PedestrianInjurySeverity", palette = injury_severity$color, popup.vars = focus_columns) +# +# tmap_save(file = "figures/dynamic_crash_maps/dynamic_crash_map.html") # generate map with leaflet ---- @@ -138,7 +138,7 @@ tag.map.title <- tags$style(HTML(" ")) title <- tags$div( - tag.map.title, HTML("Pedestrian Crashes
2017-2023") + tag.map.title, HTML("Pedestrians & Bicyclists hit by cars
2017-2023") ) tag.map.subtitle <- tags$style(HTML(" @@ -156,7 +156,7 @@ tag.map.subtitle <- tags$style(HTML(" ")) subtitle <- tags$div( - tag.map.subtitle, HTML("data from UW TOPS lab - retrieved 4/2024
per direction of the WisDOT Bureau of Transportation Safety") + tag.map.subtitle, HTML("data from UW TOPS lab - retrieved 3/2024
per direction of the WisDOT Bureau of Transportation Safety") ) leaflet() %>% @@ -193,3 +193,6 @@ leaflet() %>% # addLegendSize(position = "bottomright", color = "black", shape = "circle", values = County_Crash_Data$value.y, group = "Counties", title = "Population of County") %>% groupOptions(group ="Counties", zoomLevels = 1:9) + + +