diff --git a/howTo/openGeogAPI/apiOpenGeography-complete.R b/howTo/openGeogAPI/apiOpenGeography-complete.R
new file mode 100644
index 0000000000000000000000000000000000000000..03b056284a7ea90552742647281e22c6d2c2a07b
--- /dev/null
+++ b/howTo/openGeogAPI/apiOpenGeography-complete.R
@@ -0,0 +1,75 @@
+# get feature layer data from ONS open geography API
+
+# Really useful example of using ONS open geography portal with R
+# https://medium.com/@traffordDataLab/pushing-the-boundaries-with-the-open-geography-portal-api-4d70637bddc3
+
+library(sf)
+
+las_to_load <- c("Southampton","Portsmouth","Winchester",
+                 "Eastleigh","Isle of Wight","Fareham",
+                 "Gosport","Test Valley","East Hampshire",
+                 "Havant","New Forest","Hart","Basingstoke and Deane")
+
+# Load GHG emissions data ----
+
+library(readxl)
+
+url_to_get <- "https://assets.publishing.service.gov.uk/government/uploads/system/uploads/attachment_data/file/894787/2005-18-uk-local-regional-co2-emissions.xlsx"
+
+tempf <- tempfile(fileext = ".xlsx")
+download.file(url_to_get, tempf, method = "curl")
+
+dt <- readxl::read_xlsx(tempf, sheet = "Full dataset",skip = 1)
+
+head(dt)
+
+library(dplyr)
+dt <- dt %>%
+  filter(Name %in% las_to_load)
+
+library(reshape2)
+ghg_emissions <- melt(dt, id.vars = c("CTRY18NM/RGN18NM","Second Tier Authority","Name","Code","Year"))
+rm(dt)
+
+# Load LA geography ----
+
+# URL as API query - sometimes we don't want all boundaries
+geo_query <- "https://ons-inspire.esriuk.com/arcgis/rest/services/Administrative_Boundaries/Local_Authority_Districts_December_2018_Boundaries_UK_BGC/MapServer/0/query?where=lad18nm%20IN%20(%27Southampton%27,%27Portsmouth%27,%27Winchester%27,%27Eastleigh%27,%27Isle%20of%20Wight%27,%27Fareham%27,%27Gosport%27,%27Test%20Valley%27,%27East%20Hampshire%27,%27Havant%27,%27New%20Forest%27,%27Hart%27,%27Basingstoke%20and%20Deane%27)&outFields=lad18cd,lad18nm,long,lat&outSR=4326&f=geojson"
+
+message("Loading LA geometry from ONS Open Geography API")
+sf_data <- st_read(geo_query)
+#plot(st_geometry(sf_data))
+
+# Useful lookup spatial reference for CRS
+# https://spatialreference.org/ref/epsg/27700/
+st_coord_sys <- st_crs(sf_data) # check coord system
+st_coord_sys # current coord system EPSG: 4326 (is what leaflet wants - good)
+
+# transform the coord system if required
+if(st_coord_sys$epsg != 4326){
+  sf_data <- st_transform(sf_data, "+proj=longlat +datum=WGS84")
+}
+
+# Create map (leaflet) ----
+
+# create popup first (requires htmltools)
+library(htmltools)
+sf_data$popup_text <-
+  paste("Locial authority area code: ","<b>", sf_data$lad18cd, "</b>",
+        '<br/>', 'Local authority: ', '<b>', sf_data$lad18nm, '</b>', ' ') %>%
+  lapply(htmltools::HTML)
+
+
+library(leaflet)
+leaflet(sf_data) %>%
+  addTiles() %>%  # Add default OpenStreetMap map tiles
+  addPolygons(color = "blue", fillColor = "blue", fillOpacity = 0.2, weight = 1.5, popup = ~(lad18cd), # popups clicked
+              label = ~(popup_text),                                            # define labels
+              labelOptions = labelOptions(                                      # label options
+                style = list("font-weight" = "normal", padding = "2px 2px"),
+                direction = "auto"),
+              highlight = highlightOptions(
+                weight = 5,
+                color = "#666",
+                fillOpacity = 0.7,
+                bringToFront = TRUE))