From 986fe916afb3d4926fbcb76fe1de3287c3b72af4 Mon Sep 17 00:00:00 2001
From: Jeremy Kidwell <bots@jeremykidwell.info>
Date: Thu, 12 Oct 2023 09:14:32 +0100
Subject: [PATCH] first commit

---
 nomis.R               |  1 +
 openstreetmap_parse.R | 25 ++++++++++++++++++++++
 ordnance_survey.R     | 48 +++++++++++++++++++++++++++++++++++++++++++
 3 files changed, 74 insertions(+)
 create mode 100644 nomis.R
 create mode 100644 openstreetmap_parse.R
 create mode 100644 ordnance_survey.R

diff --git a/nomis.R b/nomis.R
new file mode 100644
index 0000000..16853c9
--- /dev/null
+++ b/nomis.R
@@ -0,0 +1 @@
+nomis.R
\ No newline at end of file
diff --git a/openstreetmap_parse.R b/openstreetmap_parse.R
new file mode 100644
index 0000000..d4ce932
--- /dev/null
+++ b/openstreetmap_parse.R
@@ -0,0 +1,25 @@
+require(sqldf) # using sqldf to filter while loading very large data sets
+require(data.table)
+require(sf) # new simplefeature data class, supercedes sp in many ways, also includes gdal functions
+
+# setup workspace
+
+setwd("/Users/kidwellj/gits/mapping_worship")
+
+if (dir.exists("data") == FALSE) {
+  dir.create("data")}
+
+# Import filtered geojson files
+# Note: preferred option here is to use R package osmdata, but the server loads generated by this extraction are too high and cause timeouts. 
+# Processing data outside R using osmium for now: https://github.com/mapcomm/cookbook/blob/master/osm_extracts/osm_extract_pow_osmium_uk.sh
+
+osm_uk_points <- st_read(system.file("pow_osm.gpkg", package = "spData"))
+
+vector_filepath = system.file("data/osm-gb-2018Aug29_pow_osm.pbf", package = "sf")
+osm_uk_points = st_read(vector_filepath)
+
+rbind(osm_uk_points, ___)
+
+# Processing handled outside R for now in QGIS
+
+# Command for identifying candidates in OSM raw csv files: grep -E 'Masjid|islamic|muslim|mosque|Madrasah|Imam' *|wc -l
\ No newline at end of file
diff --git a/ordnance_survey.R b/ordnance_survey.R
new file mode 100644
index 0000000..6886d7f
--- /dev/null
+++ b/ordnance_survey.R
@@ -0,0 +1,48 @@
+require(sf) # simplefeature data class, supercedes sp in many ways, also includes gdal functions
+
+# Data is open access, but there isn't a static link that seems to be readily 
+# accessible so I would advise that users download this directly from 
+# https://osdatahub.os.uk/downloads/open/OpenMapLocal
+# Select area: "All of Great Britain" and data format: GeoPackage
+# Fair warning, this will produce a massive file, around 3.5gb compressed and over 7gb uncompressed
+
+# This operation currently fails, but am keeping it here just for later modification
+
+if (file.exists(here("data", "opmplc_gb.gpkg", "data", "opmplc_gpkg_gb")) == FALSE) {  download.file("https://omseprd1stdstordownload.blob.core.windows.net/downloads/OpenMapLocal/2023-04/allGB/Geopackage/opmplc_gpkg_gb.zip?sv=2022-11-02&spr=https&se=2023-10-09T19%3A23%3A07Z&sr=b&sp=r&sig=83IzfszzgJxmkGGS3N7oxnpmdKvUD3JLpPQERXv0WTM%3D", destfile = here("data", "opmplc_gpkg_gb.zip"))}
+unzip("data/opmplc_gpkg_gb.zip", exdir = "data")
+
+# For academics in the UK, an alternative is to use the DigiMap service, which also provides an enhanced
+# PointX database. Sadly, this dataset is paywalled by digimap, though most academic researchers 
+# should be able to aquire it through an institutional subscription from Digimap, Ordnance Survey Data Download. 
+# The POI set is located under "Boundary and Location Data" under "Points of Interest," "select visible area" 
+# for whole UK then download.
+
+# Load in data, specifying the layer "important_building"
+# We're going to use a more sophisticated st_read command which can filter this very large file
+# h/t to https://gis.stackexchange.com/questions/341718/how-do-i-read-a-layer-from-a-gpkg-file-whilst-selecting-on-an-attribute for this approach
+
+# The feature_code for places of worship in the OpenMap product is Feature Code: 15025. This is equivalent to PointX "6340459"
+
+
+os_openmap_pow <- st_read(here("data", "opmplc_gpkg_gb", "data", "opmplc_gb.gpkg"),
+                          layer="important_building",
+                          query= "SELECT * FROM important_building WHERE feature_code = '15025';")
+
+st_write(os_openmap_pow, dsn = con, layer = "ecs",
+         overwrite = FALSE, append = FALSE)
+
+
+# These features are stored as polygons by OS so we can't export to a file as point data 
+# so CSV will end up with a blank column. FYI, if this weren't the case, here's how you'd export to CSV
+st_write(os_openmap_pow, here("data", "os_openmap_pow.csv"), layer_options = "GEOMETRY=AS_XY")
+
+# Here's how you'd export to (yuck!) shapefile:
+st_write(os_openmap_pow, here("data", "os_openmap_pow.shp"), delete_layer = TRUE)
+
+# St_wrte can also write to sql databases, and if you are working with massive complex datasets
+# it can soemtimes be more efficient to load the dataset into a PostGRES database with PostGIS extensions
+# loaded. You can then run queries from R on that database drawing in specific subsets of data using 
+# geospatial sql operations. I give an example of this in the `openstreetmapparse.R` file in this repo
+
+# And ultimately, st_write is a wrapper for GDAL, so you can export to a wide range of file types
+# including geoJSON, gpx etc etc.
\ No newline at end of file