Attempt to extract percentages of S dropout from PHE graph

library(tidyverse)

## -- Attaching packages --------------------------------------------------------------------------------------------------- tidyverse 1.3.0 --

## v ggplot2 3.3.2     v purrr   0.3.4
## v tibble  3.0.3     v dplyr   1.0.2
## v tidyr   1.1.2     v stringr 1.4.0
## v readr   1.3.1     v forcats 0.5.0

## -- Conflicts ------------------------------------------------------------------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

library(sp)
library('sf')

## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1

ltlas <-
  read_sf(
    "Local_Authority_Districts__December_2017__Boundaries_in_Great_Britain.shp"
  )

ggplot(ltlas)+geom_sf()+geom_point(aes(x=bng_e,y=bng_n))

library(png)
img <- readPNG("phe_sgtf.png")


scilly_x=104
scilly_y=843
scilly_n = 11447
scilly_e = 91327

scale_factor_x = 0.001055
scale_factor_y=-0.00124

offset_x = scilly_e - scilly_x/scale_factor_x

pos_scilly_x = (scilly_e - offset_x) * scale_factor_x

offset_y = scilly_n - scilly_y/scale_factor_y

pos_scilly = (scilly_n - offset_y) * scale_factor_y

ltlas$pixel_x = round((ltlas$bng_e - offset_x) * scale_factor_x)
ltlas$pixel_y = round( (ltlas$bng_n - offset_y) * scale_factor_y)

ltlas$color = NA
ltlas$area <- st_area(ltlas)
library(units)

## udunits system database from C:/Users/Theo/Documents/R/win-library/3.6/units/share/udunits

min_area = 36033464

units(min_area)  <- as_units("m^2")
ltlas<-ltlas %>% filter(area>min_area)
ltlas$percent =NA



deep_purple = c(62/255,16/255,85/255)
white = c(1,1,1)

color_to_percent <- function(color){
  
  full_range = deep_purple - white
  estimates = (color - white)/full_range
  if(sd(estimates)>0.07){
    return(NA)
  }
  return( mean( estimates))
  
}


for(i in 1:nrow(ltlas)){
  try_x = ltlas$pixel_x[i]
  try_y = ltlas$pixel_y[i]
  if(try_y>0){
 # print(try_x)
 #  print(try_y)
    pixel = img[try_y,try_x,]
ltlas$color[i] = list(pixel)
ltlas$percent[i] = color_to_percent(pixel)
img[try_y:(try_y+2),try_x:(try_x+2),1]=1
img[try_y:(try_y+2),try_x:(try_x+2),2]=0
img[try_y:(try_y+2),try_x:(try_x+2),3]=0
}

}







require(grDevices)
op <- par(bg = "thistle",mar = rep(0, 4))
plot(c(100, 250), c(300, 450), type = "n", xlab = "", ylab = "")

rasterImage(img, 100, 300,250, 450, interpolate = FALSE)