Cleaning images with Shiny and ImageMagick

I have a pile of old photos, literally (printed, loose or in albums and slide boxes) and figuratively (image files). Many of these are quite dirty: with dust, hairs and (particularly with scanned slides) black edges. Even fingerprint damage in places.

I also have a bunch of tools for processing photos.

But the tools just don't seem to meet my need for bulk cleaning of images. Look up guidance on how to use them, and it either seems to involve

•  using some 'heal' brush on each bit of dust (really not practical for a really dusty image), 
• or applying some 'remove dust' filter, which (for the tools I have) seems to change too much of the image, leaving non-dusty and interesting bits of the photo blurred,
• or using a modern scanner with some sort of infrared detection of dust and hairs (and re-scanning everything?).

So, time to see if I can do something that meets my need better.

Using Shiny to build a simple web application and ImageMagick to do the image processing, I think I've got quite close. The main thing I like about my approach is that you can see directly where the dust has been found, and only those pixels get changed. So you don't notice later that the faces, hair or clothes have lost focus.

The code is up at GitHub, which also has links to a shinyapp.io instance (though you're probably better running the code on your own machine).

Not perfect, but the workflow is quite quick for me now. Some of the scans done 15 or 20 years ago are jpegs, which isn't ideal. 

In future, a 'batch' version could be implemented pretty easily, but I think I need to experiment with more images before doing that.

For "tidy-select" read "column", for "data-masking" read "row"

The moment arrived for an update to himach. I thought I was solving a few issues created by newer versions of packages that himach depends on (probably due to how I'd used them in the first place). Instead, I was faced with a host of:

use of .data$ in tidyselect expressions is deprecated

warnings.

It seems that I've been using the tidyverse wrong since I first made the 2speed project that became himach, 3 years ago. I'd understood that, in a package, you need to be specific and use, say, select(.data$x) or mutate(y = .data$x + 1) so that the code can be sure you're referring to the variable x from the dataframe and not another variable. This was one of the tricky steps to get used to when moving from writing 'normal' open code, to writing a package.

But it was more subtle than that, or at least it is now, sometimes that's true and sometimes not. The tidyverse of course improves and develops, and maybe it wasn't clear to the authors then either at the time. So, updating my code wasn't just about a global search and replace, it depends on what function is being used.

There are two types of reference to variables. This seems to be a key reference.

1. <tidy-select> functions. In my code these were unnest, select, rename, across, pull. For these I did indeed have to replace .data$x with "x".  The warning message gives the right guidance. I had one case of .data[[var]]. For that the blog says use all_of(var). Instead I used {{ var  }}, because this feels right for passing a variable name as a function parameter.
2. <data-masking> functions. In my code these were mutate, group_by, filter, arrange, summarise and also (because they're within a mutate?) case_when and if_else. These cases I had to leave as .data$x. 

The logic wasn't entirely clear to me. Tidy selection is sort of manipulating columns of the dataframe, while with data-masking you're more interested in manipulating the contents. If that's the case, why does select use tidy selection, but group_by doesn't? Is it that group_by implicitly uses the contents? That has to be it.

In the end, I think it's easier to think of:

1. <column> functions, not <tidy-select>. These all manipulate the dataframe columns without caring about the rows. (Though tidyverse, I think, would like you to think of these as another sort of variable, not columns of a dataframe.)
2. <row> functions, not <data-masking>. Since these functions do different things depending on the values in each row.

Detecting objects in Images - with R, yolo and ImageMagick

I wanted to identify moving objects in video (or still) images - preferably using R as that's my tool of choice. I thought this might be challenging, but it proved to be relatively easy. 

The basic steps are:

1. Import the images from the video, discard the 'empty' ones and save the interesting ones. [R, with ImageMagick]
2. Run the interesting images through a neural net to locate and identify objects. Save the locations for some later stage of the project. [R, yolo]
3. Stitch the results back together. [R, ImageMagick]

There's an example here, applied to video from a 'camera trap' in the garden.

The rest of the blog gives the full details on how this is done.

purrr magic!

Used purrr::pmap for the first time and it's brilliant! 

Run the same function multiple times changing a bunch of parameters, and combine all the results in a single data frame. All in a single call.

Keeping code and data separate is, of course, good practice. But it can be easy to slip into mixing your code up with metadata (in this case data on where to find data and what it is). pmap makes it really easy to keep metadata separate from the code, too.

In my case, I'm pulling sets of data out of spreadsheets, where the blocks start on different rows in each sheet (don't ask!). So I just create a csv file with the parameter names on the first line, and the parameter values on the second line.

sheet,startRow,book,line

1, 53,What Might,Classics

2, 70,Runagates,Classics

3, 43,Desire,Classics

4, 45,Vocations,Classics

etc.

Read this metadata into a data frame, and then run my 'readSales' loading function, is just two lines of code:

classicSets <- read.csv(file = "data/classicSets.csv", stringsAsFactors = F)

classicSales <- pmap_dfr(classicSets, readSales)

'readSales' is just a wrapper around the excellent readxl:read_xls that cleans up my data and adds some identifying columns (book, line).

And you're done.  Thanks to Hadley Wickham and Lionel Henry!

Sustained high speeds along the Warminster Rd, Bath

In the last few posts, I was analysing number plate reader (ANPR) data from B&NES Council and BathHacked, a data activist organisation, looking at how traffic patterns might change if a clean air zone were introduced.

I've gone off on a tangent, analysing data from just two ANPRs, a mile apart on the A36 Warminster Rd, on the outskirts of Bath. They give me average speeds for some 65,000 transits over two weeks.

This is a residential 30mph zone, with a straight section, blind bend and narrow hill, though perhaps you might not think that, judging from this box plot (a smaller sample of 3 days). Each dot is a vehicle along the road. The lowest speeds are where they stop or detour (I'll put a lower limit of 5mph on to eliminate these). November 8 & 9 saw a little queueing heading into town in the morning, but these were the only cases where the speeds are not clustered around 30mph. (That's "around" - little sign of a limit in these data.)

From the data, we can split vehicles into cars, light commercial (LCVs), heavy commercial (HCVs - this is a trunk route with large lorries), public-service vehicles (PSVs) and a rare few others that we'll ignore. The scary thing is that there's no sign that trucks are travelling any more slowly than cars. Only public-service vehicles are noticeably slower, given the 3 bus stops on this section. (All days of data are included here.)

In this large survey, at peak hours there were typically 40 vehicles per hour sustaining 33mph or more over the 1 mile route. The top sustained speed was over 60mph, and one van driver managed to get 3 times into the top 20, each time doing more than 50mph.

Definitely some room for improving safety here.

Code is available on github, data from BathHacked.org.

Quick legends with chord plots

Quick legends with circlize chord plots

The circlize package from Zuguang Gu is fantastic for drawing chord plots, amongst other things. It will even allow multiple arrows between two nodes.

library(circlize)

#create quick dataframe

links <- expand.grid(c("From_A","From_B"), c("To_C","To_D", "To_E"), c("cars", "vans"))

names(links) <- c("from", "to", "type")

links$traffic <- sample.int(15, size = nrow(links))

#colours

links$colour <- ifelse(links$type == "cars", "red", "blue")

chordDiagram(links[,c("from", "to", "traffic")],

             col = links$colour,

             directional = 1, direction.type = "arrows", link.arr.type = "big.arrow")

The documentation is excellent. The catch in this case is that, if you want to add a legend, Gu refers you to yet another package ComplexHeatmap, which isn’t (currently?) available on CRAN, but via bioconductor. If, like me, you prefer not to have one package per line of code ;-) , all is not lost. There’s a quick way to get a legend, with base graphics, and another function already in circlize.

You need graphics::legend, plus circlize:add_transparency so that the colours in your legend match those used for the chords.

#for consistency

transparency <- 0.5

chordDiagram(links[,c("from", "to", "traffic")],

             col = links$colour,

             directional = 1, direction.type = "arrows", link.arr.type = "big.arrow",

             transparency = transparency)

legend(x = "bottomleft", inset = 0.05, bty = "n",

           legend = c("Cars", "Vans"),

           fill = add_transparency(c("red", "blue"), transparency),

           border = "white")

Nearly 30% of heavy commercial vehicles transiting Bath pass through Bathampton (on the A36).

We looked in earlier blogs at the overall patterns of traffic transiting Bath, using council and BathHacked data. 

In this blog we drill into the detail by type of vehicle, essentially car, and light- or heavy-commercial vehicle. We show that more transiting heavy commercial vehicles go through Bathampton on the A36 than use the Batheaston dual-carriageway.

For variety, and to show it's not very sensitive to the time threshold, in this blog we use a 45 minute threshold, and a 2-day sample (31 October and 1 November 2017). The chord diagram, which is aligned roughly geographically, looks like this.

Splitting each chord by type of vehicle, and colouring them not by starting point, but by type of vehicle gives this. Although it shows that cars are by far the biggest number of vehicles, it's hard from this to see if there are differences in the patterns.

A diagram for each could be confusing because the scales are now different: total in and out at Bathampton (Dry Arch) is near 2,000 for cars, 700 for light, and 450 for heavy commercial vehicles. 

Cars and light commercial vehicles have quite similar patterns: for light commercials there's perhaps a little higher proportion on the Swainswick-East link. But heavy commercials are very different:

• a much higher proportion are really crossing town: Bathampton - West;
• there's a much smaller proportion of traffic South-South, or South-West.

As a result, Bathampton has more than twice the share of heavy commercial transits than its share of total transits. 12.5% of all transits (3003/23841 in our sample go in or out at Bathampton (Dry Arch), compared to 28.5% (318/1116) of heavy commercials.

In fact, more transiting heavy commercial vehicles go through Bathampton on the A36 than use the Batheaston dual-carriageway.

Bathampton-Swainswick - Already more than half across the toll bridge?

As in the last two posts, I'm working through the excellent traffic data provided by BathHacked. In this post, I find the data saying that already more than half of vehicles entering Bath at Bathampton and leaving on the by-pass at Swainswick use the toll bridge. Not sure if I believe it yet, but let me know what you think.

 The chord plot showed that

1. short transits of Bath (in and out within 30 minutes) were largely split into separate East-North and South-West groupings.
2. there were about 1,800 such transits in a day past the 'Dry Arch' corner of Bathampton (where the Warminster Road becomes 40mph as you leave/enter the village)
3. The large majority of these Dry Arch transits came from or went to the bypass at Swainswick (eg to the M4).

Taking a two-day sample of ANPR data (31 Oct 17 and 1 Nov 17), the Bath data identify 130,000 individual vehicles. Here's an updated chord diagram, showing that the patterns in the list above are true for a larger 2-day sample.

The BANES Breathe project presented their modelling and plans in the Bathampton Village Hall recently. Many were dubious about the traffic projections for traffic across the toll bridge (and hence through Bathampton High Street). So naturally I want to see what the data say.

I define a 'visit' as a sequence of ANPR readings that starts or ends at the 'radial 1' ANPRs, (See the map) which include Dry Arch, and Swainswick on the by-pass.

Zooming in on visits which start or end at Dry Arch, nearly 11,000 vehicles were observed on such visits in the two days. We know from the data which direction the vehicles were heading, so we can take the traffic from the chord plot and insert more detail.

I (a) drop the 30 minute threshold (b) check where the vehicle went between Swainswick and Bathampton. There are lots of answers to (b), but we focus on whether the vehicle was seen on Warminster Road (near Trossachs Drive, so a little further into Bath) or on London Road (just beyond Cavendish Bridge)

I tried a number of ways to visualise this. A straight ggplot may be true to geography, but is hard to read! A Sankey plot is better.

But there are quite a few very small links, so it's easier to understand if we drop the links that are for fewer vehicles than 0.5% of the total. Then we get this. 

Visits start on the left and end on the right.

To give an example, of the 6,438 visits starting at Dry Arch (the orange rectangle "Bathampton_W"):

• the majority head into town passing Trossachs Drive ("Warm'rRd_W"),
• a smaller pale blue subset is next seen leaving heading North at Swainswick,
• even smaller subsets leave again at Dry Arch ("Bathampton_E"): commuters, shoppers or day-trippers,
• or are picked up at London Road (maybe the ANPR missed them at Trossachs).

The sharp-eyed reader will have noticed more visits coming into WarminsterRd_W than leaving: the omissions are those that go on to the car parks (in which case they may account for more visits leaving WarminsterRd_E than arrive), or leave by other routes, so aren't of interest for the current analysis of Bathampton-Swainswick.

We could get more detail, but we are interested in the pale blues: the hypothesis being that these go over the toll bridge. Of the 6,438 visits starting at Dry Arch, 1,126 (17.5%) might have gone over the toll bridge. This is an upper estimate, which could more precisely be put that we know that 82.5% did not.

There are 996 possible toll bridge crossings in the other direction, making for around 1,000 per day, in total in the two directions. 

This 17.5%, however, is not the percentage we are looking for. Instead we need a more interesting one. Note the share of dark and light blue at Swainswick_N. Of our 2,077 total visits coming in at Dry Arch and leaving at Swainswick, 1,126 have possibly crossed the bridge. So 54% are already crossing the toll bridge.

This is a surprisingly high percentage to me, given 

1. that most lorries cannot use the bridge
2. a toll is a toll
3. the queues at the toll bridge are unpredictable
4. we only had about 1,300 30-minute transits (see chord diagram), so what happens to the other 2077-1300 ~ 750?

So this definitely needs further checking. The fact that the situation is similar in these data in the opposite direction gives some confirmation, but further validation, for example looking at travel times and vehicle types is for another day.

Notes

Implicitly, this analysis excludes

• around 5% where the ANPR did not identify the vehicle (I don't know the technology, but perhaps the number plate was obscured by another vehicle, or was in an unknown format), 
• where only a single ANPR reading was made in the 2 days.

and hat-tip to the networkD3 package authors, and all the other contributors to R

Bathampton - transiting traffic

Building on the previous blog post, I've extracted some 13,000 vehicle transits of Bath on 31 October 2017 from data provided by BathHacked. (These are all the entries and exit again within 30 minutes, through the 'radial 1' on the map shown before.)

A chord plot, using the great circlize package in R is a good way to show these flows. The counts around the edge are transit entries or exits - so the total is 26,000. Direction is indicated by the slight point/arrowhead.

You can see a big influence of the Batheaston bypass in the flow between North (Lansdown and Swainswick (dual carriageway)) and East (Box Road & Bathford). That's expected.

But I'm really surprised that the very large majority of vehicles entering or leaving Bathampton (near Dry Arch) and transiting Bath are also heading to, or coming from the North. Around 700 heading North, and 1,000 heading South. Very little traffic transiting and heading West, for example.

Can this be correct? More checking needed, I think.

(Note that one vehicle might generate more than one transit in the day

How long will I taxi at at airport?

CODA Taxi Times

Eurocontrol/CODA has been publishing aircraft taxi times by airport for some years. These are calculated, to the nearest minute, from flight-by-flight data provided to CODA by airlines and airports. Airports are included where CODA receives data on more than 100 flights, so it covers airports small and large, mostly in Europe, but some non-European ones with direct flights to Europe.

Taxi out is the time from pushback from the gate, to take off. Taxi in the time from landing to on-blocks at the gate. (With some minor variation for remote stands..)

I've now published in github some R code for scraping all of these data into a single, tidy dataset. The graph gives an example for 5 big European airports: showing how taxi-in times are shorter than taxi-out, but that there has been little variation over the years (or for that matter between seasons).

After each IATA season (Summer = end March-end October, Winter), CODA publishes 4 reports: taxi-in times, taxi-out times, and the same, but split by ICAO wake turbulence category (WTC). The files were pdf for the first few years, now xlsx.

This project was for me a learning exercise in web-scraping, and in using git and github. 

However, the pdf reports by WTC were difficult to read because of the row alignment in the tables, so these are not currently included. Perhaps tabulizer would handle this, but I gave up after wasting some hours trying to get its java to work on my machine, and parsed the tables manually from the text, instead!

SWDChallenge - Dot Plot

Eurocontrol's CODA team publishes regular stats on airline punctuality in Europe, based on data provided directly by airlines.

Some of the graphics are showing their age, so the dot plot challenge was a good excuse to try something slightly different.

I picked the 'punctuality' graph. CODA doesn't name individual airlines, because this is about performance benchmarking, not pointing the finger. Airlines are told their own punctuality, then the graph allows them to compare against the other largest airlines. 

I wanted to achieve two things:

1) Make it easier to compare quarters (currently an airline will be in one horizontal position in one quarter, and another in the other).

2) Highlight the airlines with the most flights. Each airline's punctuality challenges are different, depending on where they fly, their fleet, whether they offer connecting flights etc. Improving punctuality provides a better service to passengers, whatever the size of airline - but larger airlines have more effect on the average performance of the whole network.

The data look something like this (with airline names anonymised..), after some ranks have been added, with a mutate_at(vars(Flt1,Punc1), funs(rank = min_rank(desc(.)))).  I fiddled with the data a bit to get some good tests (of decline, improvement, and stable), but in any case they are not final or official.

I tried two options. 

In both, I went for using triangles and tails to indicate the change since the previous quarter. But the triangles tend to point to the wrong place. Rather than work out how to nudge the triangles, I used a dot to emphasise the actual value.

The first option used size to distinguish larger airlines from smaller. I tried various options for dot, triangle size and transparency (alpha), but in the end, felt that there was just too much going on, and the big airlines did not stand out.

So I switched to using colour instead. It's still quite busy, but I think the result is better, and a step forward from the current 'legacy' graphs.

For the record, this was the code.

#calculate some values for the plot

w <- punc %>% 

  #edit some values to get good test

  mutate(Punc0 = if_else(Punc1_rank %in% c(14, 32,44), Punc1, Punc0)) %>% 

  #calculate shapes and size

  mutate(col = if_else(Flt1_rank <= bigRank, "black", "cyan3"),

         shape = case_when(

           Punc1 == Punc0 ~ 1, #circle

           Punc1 >= Punc0 ~ 24, #up arrow

           TRUE           ~ 25))

#where to label

annoX <- min(w$Punc1)

ggplot(w) + 

  #last quarter

  geom_segment(aes(x = Punc1_rank, y=Punc0, xend= Punc1_rank, yend=Punc1, colour = col), alpha=0.5)+ 

  #this quarter - direction triangles

  geom_point(aes(x = Punc1_rank, y=Punc1, shape = shape, colour = col), size = 3) +

  #this quarter - points

  geom_point(aes(x = Punc1_rank, y=Punc1, colour = col), size = 0.7) +

  annotate("text", x = 1, y = annoX, size = 2.5, 

           hjust = 0, vjust = 0,

           label = paste("Black = Top 10 airlines by flights in",q1,

                          "\nArrow tail = punctuality in",q0)) +

  scale_shape_identity() + #use the value directly as a shape

  scale_colour_identity() + #use the value directly as a size

  labs(x="Rank of airline by punctuality (best to worst). Top 50 airlines by flights are shown.",

       y=paste("Arrival Punctuality in",q1,"(Delay < 15 minutes)")) +

  scale_x_continuous(breaks = c(1, seq(5, maxRank, by = 5)),

                     minor_breaks = NULL) +

  scale_y_continuous(labels= scales::percent) +

  theme_minimal()

Web scraping for cycling results - Here be pirates!

In the previous blog post, I went through the steps of downloading data from a number of sequential pages, and saving them on your local system as separate html files. 

I use the tidyverse, and also package chron, since I have data which are pure elapsed times.

I then used the rvest package to parse these files. You need to inspect the html to see the structure of the data as now saved in the .html file. You can just open the file in RStudio, for example, and scroll down until you see some of the data.

In the case of this cycling data, there are rows 'tr' alternating of class 'Odd' and 'Even'. Then within each row, there are a number of data elements, 'td'.

We mimic this structure with two calls to 'html_nodes()', the first picking out the rows by giving the class (with an extra . in front), the second pulling out the data elements. Then we use html_text to strip off the extraneous text.

 
library(tidyverse)
library(rvest) #for web-scraping
library(chron) #for times

allRes <- unlist(lapply(0:4, function(x){
  # read in the html
  res <- read_html(paste0("race",x,".html")) 
  # pick the data we want - this comes with the column headings
  z <- unlist(lapply(c(".Odd",".Even"), 
                     function(x){res %>% 
                         html_nodes(x) %>% 
                         html_nodes("td") %>% 
                         html_text(trim = TRUE)}))                   
}))
  

This is partly a hack in that I'm assuming I know what the files are, and that they have indices 0 to 4. But that's true here, if not elegant.

Finally, format all these data into a dataframe, via a matrix, for clarity. Everything arrives as strings, so I name and change the formats to something more meaningful.

#convert to dataframe
z <- matrix(allRes, byrow=TRUE, ncol=14)
colnames(z) <- c("Pos","ID","Sex","Name","Age","Nat","Team","Time","Tkm","AveSpd","CatPos","Cat","City","x")
results <- as_tibble(z) %>% 
  #reformat some columns
  mutate_at(.vars = vars(Pos, ID, Age, CatPos),
            .funs = funs(as.integer)
         ) %>% 
  #remove non-finishers
  filter(!is.na(Pos)) %>% 
  mutate(AveSpd = as.numeric(AveSpd),
         # make time always HMS
         Time = times(if_else(nchar(Time)<6, paste0("0:",Time), Time)),
         #complete sex
         Sex = if_else(Sex=="F",Sex,"M")) %>% 
  #we still get some NAs, which are the ones whose age is unknown, so their CatPos is also unknown
  drop_na() 
  # manual ordering
  arrange(Cat, AveSpd)

After that, it's hard to resist having a look at the data. I was keen to try out the pirate plot for the first time, from package yarrr. This is an easy way to view data in a number of categories: viewing all of the data points (as jittered points), a distribution density for each category, and some descriptive statistics. For the last of these, I chose the option 'iqr' for interquartile range, which gives me something looking like a box plot on top of the data, especially when coupled with the median function for the central line.

Here, the categories are competition age groups.

library(yarrr) #for pirate plot
#plot of speed
pp <- pirateplot(formula = AveSpd ~ Cat ,
    data = results,
    avg.line.fun = median,
    ylab = 'Average Speed (kmh)',
    xlab = "Competition Group",
    inf.method = 'iqr',
    theme = 2,
    point.o = 0.6
    ) 

The resulting plot is lovely: