latlon2poa2016

#' @return Character vector of postcodes on 2016 basis
latlon2poa <- function(lat, lon, shapefile.dir = Sys.getenv("POA_2016_SHAPEFILE_DIR")) {
  stopifnot(!anyNA(lat), !anyNA(lon), is.numeric(lat), is.numeric(lon), 
            # Basic check for lat=lon mixup
            max(lat) <= 90)
  if (!dir.exists(shapefile.dir)) {
    tempf <- tempfile()
    dir.create(tempf)
    exts <- c("dbf", "prj", "shp", "shx", "xml")
    # Base URL
    ub <- "https://github.com/HughParsonage/POA2016/raw/master/shp/1270055003_poa_2016_aust_shape/POA_2016_AUST."
    for (ext in exts) {
      s <- download.file(url = paste0(ub, ext),
                         mode = "wb", 
                         destfile = file.path(tempf, paste0("POA_2016_AUST.", ext)), 
                         quiet = TRUE)
      if (s) {
        stop("Error #", s, " downloading\n\t", paste0(ub, ext))
      }
    }
    
    shapefile.dir <- tempf
    # to avoid downloading multiply same session
    Sys.setenv("POA_2016_SHAPEFILE_DIR" = shapefile.dir)
  }
  # rgdal expects no trailing slash
  if (endsWith(shapefile.dir, "/")) {
    shapefile.dir <- sub("/$", "", shapefile.dir)
  }
  suppressWarnings({
    # Warning messages:
    #   1: In OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS = dumpSRS,  :
    #  Discarded datum Geocentric_Datum_of_Australia_1994 in CRS definition: +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs
    #                       2: In rgdal::readOGR(shapefile.dir) : Dropping null geometries: 2669, 2670
    shapefile <- rgdal::readOGR(shapefile.dir, verbose = FALSE)
  })
  
  points <- sp::SpatialPoints(coords = sp::coordinates(data.frame(x = lon, y = lat)),
                              proj4string = shapefile@proj4string)
  out <- sp::over(points, shapefile)
  out$POA_NAME16
}