From 2c4ecf1757774a3512ca2cccfe6376341d227f7b Mon Sep 17 00:00:00 2001 From: Amy Mikhail Date: Thu, 29 Jun 2023 18:06:06 +0100 Subject: [PATCH] Added osmdata rooftop id to Population Estimates --- Population Estimates/data/random_points.gpx | 154 ------- Population Estimates/gis_sampling.Rmd | 381 ++++++++++++----- Population Estimates/gis_sampling.html | 428 +++++++++++++------- docs/gis_sampling.html | 428 +++++++++++++------- 4 files changed, 860 insertions(+), 531 deletions(-) delete mode 100644 Population Estimates/data/random_points.gpx diff --git a/Population Estimates/data/random_points.gpx b/Population Estimates/data/random_points.gpx deleted file mode 100644 index 120e928..0000000 --- a/Population Estimates/data/random_points.gpx +++ /dev/null @@ -1,154 +0,0 @@ - - - - - point-1 - - - point-2 - - - point-3 - - - point-4 - - - point-5 - - - point-6 - - - point-7 - - - point-8 - - - point-9 - - - point-10 - - - point-11 - - - point-12 - - - point-13 - - - point-14 - - - point-15 - - - point-16 - - - point-17 - - - point-18 - - - point-19 - - - point-20 - - - point-21 - - - point-22 - - - point-23 - - - point-24 - - - point-25 - - - point-26 - - - point-27 - - - point-28 - - - point-29 - - - point-30 - - - point-31 - - - point-32 - - - point-33 - - - point-34 - - - point-35 - - - point-36 - - - point-37 - - - point-38 - - - point-39 - - - point-40 - - - point-41 - - - point-42 - - - point-43 - - - point-44 - - - point-45 - - - point-46 - - - point-47 - - - point-48 - - - point-49 - - - point-50 - - diff --git a/Population Estimates/gis_sampling.Rmd b/Population Estimates/gis_sampling.Rmd index 07d7ad4..103b229 100644 --- a/Population Estimates/gis_sampling.Rmd +++ b/Population Estimates/gis_sampling.Rmd @@ -1,7 +1,7 @@ --- title: "Spatial sampling and population estimates with R" author: "Alexander Spina and Alexandre Blake (Applied Epi)" -date: "July 2022" +date: "July 2023" output: html_document: code_folding: show @@ -62,6 +62,7 @@ pacman::p_load( leaflet, # making interactive maps mapedit, # interactively drawing points and points on maps sf, # working with spatial data + osmdata, # getting data from OpenStreetMap here, # file paths relative to R project root folder dplyr, # data management rio # import data @@ -173,16 +174,19 @@ in to the leaflet map below with your mouse. ```{r camp_location} kario_camp <- leaflet() %>% - # start an interactive plot (of the whole world) - addTiles() %>% - # center the leaflet map on the middle of our area of interest - setView( - lng = 26.208, - lat = 11.1496, - zoom = 15) %>% - # add in background images of the area from open street map - # we use the humanitarian version for this setting as there is more info - addProviderTiles("OpenStreetMap.HOT") + + # start an interactive plot (of the whole world) + addTiles() %>% + + # center the leaflet map on the middle of our area of interest + setView(lng = 26.208, + lat = 11.1496, + zoom = 15) %>% + + # Add in background images of the area from open street map + # We use the humanitarian version for this setting as there is more info + addProviderTiles("OpenStreetMap.HOT") + ``` You can view and interact with your map simply by calling the object you @@ -190,7 +194,7 @@ saved to, as below. ```{r camp_view} -# view your map interactively in browser +# View your map interactively in browser kario_camp ``` @@ -210,14 +214,14 @@ need additional help. ```{r camp_polygon_draw, eval = FALSE} -# interactively draw a shape on your map -# your polygon is saved in the object "polygon" (which is a list) -# more specifically in the slot "finished" of the object polygon. -polygon <- editMap(kario_camp) +# Interactively draw a shape on your map +# Your polygon is saved in the object "polygon" (which is a list) +# More specifically in the slot "finished" of the object polygon. +kario_polygon <- editMap(kario_camp) -# access your polygon -polygon <- polygon$finished +# Access your polygon +kario_polygon <- kario_polygon$finished ``` #### Task C: Read in an existing polygon (optional) @@ -239,9 +243,9 @@ knows to pull in the others. ```{r camp_polygon_read} # read in your shapefile -polygon <- read_sf(here::here("data", - "SDN_Kario_17Feb2019_polygon", - "SDN_Kario_17Feb2019_polygon.shp")) +kario_polygon <- read_sf(here::here("data", + "SDN_Kario_17Feb2019_polygon", + "SDN_Kario_17Feb2019_polygon.shp")) ``` @@ -250,16 +254,20 @@ You could then view the polygon to your existing map. ```{r camp_polygon_plot} kario_camp <- kario_camp %>% - # plot your polygon on top of the existing tiles + + # Plot your polygon on top of the existing tiles addPolygons( - # specify your sf object - data = polygon, - # color the border red + + # Specify your sf object + data = kario_polygon, + + # Color the border red color = "red", - # make the middle see-through + + # Make the middle see-through fillOpacity = 0) -# view your map +# View the map: kario_camp ``` @@ -294,24 +302,133 @@ need additional help. ```{r household_points, eval = FALSE} -# interactively draw a shape on your map -# your points are saved in the object "households" (which is a list) -# more specifically in the slot "finished" of the object households -households <- editMap(kario_camp) +# Interactively draw a shape on your map +# Your points are saved in the object "households" (which is a list) +# More specifically in the slot "finished" of the object households +households <- editMap(kario_map) -## access your points +# Access your points households <- households$finished ``` -#### Task B: Read in existing rooftop points (optional) + +#### Task B: Automatically identify rooftops (optional) *This is an alternative step that replaces marking roofs in the previous task* -The full roofs identification of the refugee camp was in fact carried -out semi-automatically, manually validated and saved in a shapefile. As -with the polygon above we can use the `sf` package to read in the points -to R. +It is possible to identify each rooftop in the camp automatically, using an OpenStreetMap feature that identifies buildings from satellite imagery. We can do this in R with the package `osmdata`. + +First, we need to define a boundary box within which all the building data will be retrieved. We can get the boundary box for the Kario camp polygon as below: + +```{r kario_boundarybox} + +# Get a boundary box for Kario camp from the shapefile: +kario_bbox <- sf::st_bbox(kario_polygon) + +``` + + +Now we can retrieve the geometries (polygons and points) for all the buildings that fall within this boundary box. We use the `opq` function to query OpenStreetMap and select "building" as the feature that we want to retrieve: + +```{r kario_buildings} + +# Get the builing points within the boundary box: +buildings <- osmdata::opq(bbox = kario_bbox) %>% + + # Get the buildings in this boundary box: + osmdata::add_osm_feature(key = "building") %>% + + # Convert the building points to an osmdata_sf object: + osmdata::osmdata_sf() + + +``` + + +If you inspect the contents of `buildings` you can see that it is a list of different objects, including two data.frames: `osm_points` and `osm_polygons`. The coordinates in `osm_points` define the polygons for each building. Note that this means there are multiple points per building (one point for each corner in the shape that makes up the building). We can already see the outlines of the buildings on the map tiles however, so what we really need are a single pair of coordinates for each building (to represent one household). We can get this by calculating the centroid for each building polygon: + +```{r building_centroids} + +# Select the osm_polygons data.frame: +households_osm <- buildings$osm_polygons %>% + + # Convert it to a simple features object: + sf::st_as_sf() %>% + + # Get centroids (points) for each building polygon: + mutate(points = sf::st_centroid(geometry)) + + +``` + + +Now we can add the points representing each household to the map: + +```{r map_buildings} + +# Update the map: +kario_map <- kario_camp %>% + + # Add the points representing buildings to the map: + addCircleMarkers( + + # Get building points to plot from the points column: + data = households_osm$points, + + # Specify the radius of the circle markers: + radius = 5, + + ) + + +# View the updated map: +kario_map + +``` + + +Currently some of the buiding points are located outside the Kario camp polygon - this is because the query we used to fetch the building data used the boundary box of the polygon, not the polygon itself. We can remove the points outside the polygon by trimming the data with the `st_intersection()` function. This takes two arguments - x (the points data that we want to trim) and y (the polygon data that we want to use for the trimming). + + +```{r remove_excess_buildings} + +households_osm <- households_osm %>% + + # Trim to points which are inside the Kario camp polygon: + sf::st_intersection(y = kario_polygon) %>% + + # St_intersection adds an extra id column we don't need: + select(-id) + + +``` + + +Now we can redraw the map and check it to make sure the extra points have been removed: + +```{r kario_map_noexcess} + +# Update the map of Kario camp: +kario_map <- kario_camp %>% + + # Plot the points representing buildings ontop: + addCircleMarkers( + data = households_osm$points, + radius = 5, + ) + + +# View the map: +kario_map + +``` + +There is still a problem, however. The automated process has identified each building from satellite imagery, but not all of these buildings are occupied, and some of them are not households but have other functions (such as central storage facilities or bathroom blocks). In addition, some new families have arrived since the last satellite image was taken. + +To ensure that each point represents an occupied household, the automatically generated points were validated by a field survey team and saved in a shapefile. + +As with the polygon above we can use the `sf` package to read in the points to R. An additional complication here is that the points are not stored as longitude and latitude, but as x and y coordinates. For this we use the @@ -322,26 +439,44 @@ Maps). ```{r camp_points_read} -# read in your shapefile -households <- read_sf(here::here("data", - "SDN_Kario_17Feb2019_roofs", - "SDN_Kario_17Feb2019_roofs.shp")) %>% - # re-project x and y coordinates to be longitude and latitude. - # note that crs 4326 is a code for WGS84 +# Read in your shapefile +households_validated <- sf::read_sf( + here::here("data", + "SDN_Kario_17Feb2019_roofs", + "SDN_Kario_17Feb2019_roofs.shp")) %>% + + # Re-project x and y coordinates to be longitude and latitude. + # Note that crs 4326 is a code for WGS84 st_transform(crs = 4326) + ``` -You could then view the polygon to your existing map. + +Once again, we can redraw the map, this time using the manually validated points from the shapefile: + ```{r camp_points_plot} -kario_camp %>% - # add circles for each point - # (zoom in to see how they fit with rooftops) - addCircleMarkers(data = households) +# Update the map of Kario camp: +kario_map <- kario_camp %>% + + # Plot the points representing buildings ontop: + addCircleMarkers( + data = households_validated, + radius = 5, + ) + + +# View the map: +kario_map + ``` + + + + ## Step 3: *Sample* ### Objective @@ -365,7 +500,7 @@ Kario survey, 400 roofs were sampled in order to do Roof-Simple-Random. ```{r sample_points} # select 50 random points from your polygon (i.e. not necessarily households) -random_points <- st_sample(polygon, 50) %>% +random_points <- st_sample(kario_polygon, 50) %>% # change back to a simple features dataframe (sampling made it a "collection") st_as_sf() @@ -374,7 +509,7 @@ random_points <- st_sample(polygon, 50) %>% kario_camp %>% # add circles for each point # (zoom in to see how they fit with rooftops) - addCircleMarkers(data = random_points, color = "green") + addCircleMarkers(data = random_points, color = "green", radius = 5) ``` @@ -390,35 +525,64 @@ does not account for this. While it is not the case in this tutorial, in other situations you might have more than one polygon and need to sample roofs separately for each polygon. -Notice here when you visualise that the points are much closer to where -the population is actually situated in the camp (because you are -sampling households rather than any space within the camp boundary) - -this likely makes conducting your survey logistically simpler. In -addition, you only need to survey half as many houses compared to -T-Square (because you don't need to find the next closest house). +We will first try to sample from the OSM household (unvalidated) data: + +```{r sample_roofs_osm} + +# Select 50 random roofs from the household_osm data: +random_roofs_osm <- st_sample(households_osm, 50) + +# Visualise the selected points on a map: +kario_camp %>% + + # Add circles for each point + addCircleMarkers(data = random_roofs_osm, + color = "orange", + radius = 5) -```{r sample_roofs} -# select 50 random points from your households -random_roofs <- st_sample(households, 50) %>% - # transform back to points - # (this comes out as a multipoint object, which is too complex) +``` + + +For comparison, we can also select 50 random points from the validated data. Note that this data requires some extra formatting before we can visualise it: + +```{r sample_roofs_validated} +# Select 50 random points from your households +random_roofs_val <- st_sample(households_validated, 50) %>% + + # Transform back to points + # (this comes out as a multi-point object, which is too complex) st_cast("POINT") %>% - # set the coordinate reference system for the points - # note that crs 4326 is a code for WGS84 + + # Set the coordinate reference system for the points + # Note that crs 4326 is a code for WGS84 st_set_crs(4326) %>% - # change back to a simple features dataframe (multipoint made it a "collection") + + # Change back to a simple features data.frame (multi-point made it a "collection") st_as_sf() # Visualise the selected points kario_camp %>% - # add circles for each point - # (zoom in to see how they fit with rooftops) - addCircleMarkers(data = random_roofs, color = "orange") + + # Add circles for each point + addCircleMarkers(data = random_roofs_val, + color = "orange", + radius = 5) + + ``` + +Notice here when you visualise that the points are much closer to where +the population is actually situated in the camp (because you are +sampling households rather than any space within the camp boundary) - +this likely makes conducting your survey logistically simpler. In +addition, you only need to survey half as many houses compared to +T-Square (because you don't need to find the next closest house). + + ## Step 4: *Export* ### Objective @@ -437,19 +601,25 @@ either. To save our points we use the `st_write()` function from the ```{r export_points_gpx, results = 'hide'} -# export sampled points to gpx +# Export sampled points to gpx: random_points %>% + # gpx has a specific format which it accepts - # to meet this format we have to give a name to label each of the points - # note that this needs to be a lower case "name" - # we consecutively number the points by combining "point" with the number of the row + # To meet this format we have to give a name to label each of the points + # Note that this needs to be a lower case "name" + # We consecutively number the points by combining "point" with the number of the row mutate(name = paste0("point-", 1:n())) %>% + + # Now we can write the labelled points to a .gpx file: st_write( - # define the file path and extension (gpx) - dsn = here::here("data", "random_points.gpx"), - # define which format to use (this could be guessed by the function) + + # Define the file path and extension (gpx) + dsn = here::here("data", "random_points.gpx"), + + # Define which format to use (this could be guessed by the function) driver = "GPX", - # overwrite any previously existing files. + + # Overwrite any previously existing files. delete_dsn = TRUE) @@ -460,35 +630,47 @@ We could similarly use `st_write()` to save to a `.shp` or a `.kml` ```{r export_points_other, eval = FALSE} -# export sampled points to shp +# Export sampled points to shp random_points %>% - # it can be useful to give a name to label each of the points - # we consecutively number the points by combining "point" with the number of the row + + # It can still be useful to give a name to label each of the points + # We consecutively number the points by combining "point" with the number of the row mutate(Name = paste0("point-", 1:n())) %>% - st_write( - # define the file path and extension (kml) - dsn = here::here("data", "random_points.shp"), - # overwrite any previously existing files. + + # Now we can write the labelled points to a shapefile: + st_write( + + # Define the file path and extension (kml) + dsn = here::here("data", "random_points.shp"), + + # Overwrite any previously existing files. delete_dsn = TRUE) -# export sampled points to kml +# Export sampled points to kml random_points %>% + # kml has a specific format which it accepts - # to meet this format we have to give a name to label each of the points - # note that this needs to be a lower case "Name" - # we consecutively number the points combining "point" with the number of the row + # To meet this format we have to give a name to label each of the points + # Note that this needs to be a lower case "Name" + # We consecutively number the points combining "point" with the number of the row mutate(Name = paste0("point-", 1:n())) %>% - st_write( - # define the file path and extension (kml) - dsn = here::here("data", "random_points.kml"), - # define which format to use (this could be guessed by the function) + + # Now we can write the labelled points to a kml file: + st_write( + + # Define the file path and extension (kml) + dsn = here::here("data", "random_points.kml"), + + # Define which format to use (this could be guessed by the function) driver = "kml", - # overwrite any previously existing files. + + # Overwrite any previously existing files. delete_dsn = TRUE) ``` + #### Task B: Transfer the sampled points 1. To a mobile phone via USB cable @@ -518,12 +700,13 @@ and open the Garmin eTrex device until the .../GPX/... folder.\ Copy-paste the saved gpx file from your computer hard drive to the GPX folder of the Garmin eTrex device. -## Step 4: *Calculate* + +## Step 5: *Calculate* ### Objective The objective of this step is to calculate your population estimates -using T-square and then Roof-Simple-Random. +using T-square and then the Roof-Simple-Random method. #### Task A: Enter T-Square data @@ -599,7 +782,7 @@ our case the area of Kario camp (in *m^2^*). We can do this using the ## calculate the area of within the polygon and return as a numeric ## otherwise it is returned as a measurement (with m2 on the end) -camp_area <- st_area(polygon) %>% +camp_area <- st_area(kario_polygon) %>% as.numeric() ## view your numbers @@ -643,7 +826,7 @@ separate point objects. first_point <- survey_data %>% select(gps_point) %>% ## split in to a latitude and longitude column - tidyr::separate(gps_point,into = c("lat","lon"), sep=", ") %>% + tidyr::separate(gps_point,into = c("lat","lon"), sep = ", ") %>% ## change columns to numeric mutate(across(c(lat, lon), as.numeric)) %>% ## set as sf points and define the crs @@ -652,14 +835,14 @@ first_point <- survey_data %>% ## extract the first randomly sampled points second_point <- survey_data %>% select(nearest_house_gps_point) %>% - tidyr::separate(nearest_house_gps_point,into = c("lat","lon"), sep=", ") %>% + tidyr::separate(nearest_house_gps_point,into = c("lat","lon"), sep = ", ") %>% mutate(across(c(lat, lon), as.numeric)) %>% st_as_sf(coords = c("lat","lon"), crs = 4326) ## extract the first randomly sampled points third_point <- survey_data %>% select(second_house_gps_point) %>% - tidyr::separate(second_house_gps_point,into = c("lat","lon"), sep=", ") %>% + tidyr::separate(second_house_gps_point,into = c("lat","lon"), sep = ", ") %>% mutate(across(c(lat, lon), as.numeric)) %>% st_as_sf(coords = c("lat","lon"), crs = 4326) @@ -803,7 +986,7 @@ averages <- survey_data %>% summarise(across(c(under_five, five_plus, overall), mean)) ## pull the number of households from the shapefile -num_hh <- nrow(households) +num_hh <- nrow(households_validated) ## multiply the averages by the number of houses population <- averages * num_hh @@ -835,7 +1018,7 @@ You can read in these files with the following code. ```{r read_gpx, eval = FALSE} -## read in a track from osmand +## Read in a track from osmAnd polygon <- st_read( dsn = here::here("data", "track.gpx") ) diff --git a/Population Estimates/gis_sampling.html b/Population Estimates/gis_sampling.html index 082e72d..7bab766 100644 --- a/Population Estimates/gis_sampling.html +++ b/Population Estimates/gis_sampling.html @@ -1667,20 +1667,18 @@ document.body.style.height = "100%"; document.documentElement.style.width = "100%"; document.documentElement.style.height = "100%"; - if (cel) { - cel.style.position = "absolute"; - var pad = unpackPadding(sizing.padding); - cel.style.top = pad.top + "px"; - cel.style.right = pad.right + "px"; - cel.style.bottom = pad.bottom + "px"; - cel.style.left = pad.left + "px"; - el.style.width = "100%"; - el.style.height = "100%"; - } + cel.style.position = "absolute"; + var pad = unpackPadding(sizing.padding); + cel.style.top = pad.top + "px"; + cel.style.right = pad.right + "px"; + cel.style.bottom = pad.bottom + "px"; + cel.style.left = pad.left + "px"; + el.style.width = "100%"; + el.style.height = "100%"; return { - getWidth: function() { return cel.offsetWidth; }, - getHeight: function() { return cel.offsetHeight; } + getWidth: function() { return cel.getBoundingClientRect().width; }, + getHeight: function() { return cel.getBoundingClientRect().height; } }; } else { @@ -1688,8 +1686,8 @@ el.style.height = px(sizing.height); return { - getWidth: function() { return el.offsetWidth; }, - getHeight: function() { return el.offsetHeight; } + getWidth: function() { return cel.getBoundingClientRect().width; }, + getHeight: function() { return cel.getBoundingClientRect().height; } }; } } @@ -1913,8 +1911,8 @@ elementData(el, "initialized", true); if (bindingDef.initialize) { - var result = bindingDef.initialize(el, el.offsetWidth, - el.offsetHeight); + var rect = el.getBoundingClientRect(); + var result = bindingDef.initialize(el, rect.width, rect.height); elementData(el, "init_result", result); } } @@ -1956,29 +1954,30 @@ forEach(matches, function(el) { var sizeObj = initSizing(el, binding); + var getSize = function(el) { + if (sizeObj) { + return {w: sizeObj.getWidth(), h: sizeObj.getHeight()} + } else { + var rect = el.getBoundingClientRect(); + return {w: rect.width, h: rect.height} + } + }; + if (hasClass(el, "html-widget-static-bound")) return; el.className = el.className + " html-widget-static-bound"; var initResult; if (binding.initialize) { - initResult = binding.initialize(el, - sizeObj ? sizeObj.getWidth() : el.offsetWidth, - sizeObj ? sizeObj.getHeight() : el.offsetHeight - ); + var size = getSize(el); + initResult = binding.initialize(el, size.w, size.h); elementData(el, "init_result", initResult); } if (binding.resize) { - var lastSize = { - w: sizeObj ? sizeObj.getWidth() : el.offsetWidth, - h: sizeObj ? sizeObj.getHeight() : el.offsetHeight - }; + var lastSize = getSize(el); var resizeHandler = function(e) { - var size = { - w: sizeObj ? sizeObj.getWidth() : el.offsetWidth, - h: sizeObj ? sizeObj.getHeight() : el.offsetHeight - }; + var size = getSize(el); if (size.w === 0 && size.h === 0) return; if (size.w === lastSize.w && size.h === lastSize.h) @@ -2280,7 +2279,6 @@ return result; } })(); - +
+

Task B: Draw the polygon

@@ -6655,14 +6656,14 @@

Task B: Draw the polygon

done, click on “Finish”. Then click on “Done” on the bottom right. Have a quick look at the mapedit website if you need additional help.

-
# interactively draw a shape on your map
-# your polygon is saved in the object "polygon" (which is a list)
-# more specifically in the slot "finished" of the object polygon.
-polygon <- editMap(kario_camp)
+
# Interactively draw a shape on your map
+# Your polygon is saved in the object "polygon" (which is a list)
+# More specifically in the slot "finished" of the object polygon.
+kario_polygon <- editMap(kario_camp)
 
 
-# access your polygon 
-polygon <- polygon$finished

+# Access your polygon 
+kario_polygon <- kario_polygon$finished

Task C: Read in an existing polygon (optional)

@@ -6678,24 +6679,28 @@

Task C: Read in an existing polygon (optional)

a shapefile consists of many more files than just the .shp, if we tell R where that is, it then knows to pull in the others.

# read in your shapefile
-polygon <- read_sf(here::here("data", 
-                              "SDN_Kario_17Feb2019_polygon", 
-                              "SDN_Kario_17Feb2019_polygon.shp"))
+kario_polygon <- read_sf(here::here("data", + "SDN_Kario_17Feb2019_polygon", + "SDN_Kario_17Feb2019_polygon.shp"))

You could then view the polygon to your existing map.

kario_camp <- kario_camp %>% 
-  # plot your polygon on top of the existing tiles
+  
+  # Plot your polygon on top of the existing tiles
   addPolygons(
-    # specify your sf object
-        data  = polygon,
-        # color the border red
+    
+    # Specify your sf object
+        data  = kario_polygon,
+        
+        # Color the border red
         color = "red",
-        # make the middle see-through
+        
+        # Make the middle see-through
         fillOpacity = 0)
 
-# view your map
+# View the map:
 kario_camp
-
- +
+
@@ -6723,42 +6728,148 @@

Task A: Mark the roofs

“Finish”. Then click on “Done” on the bottom right. Have a quick look at the mapedit website if you need additional help.

-
# interactively draw a shape on your map
-# your points are saved in the object "households" (which is a list)
-# more specifically in the slot "finished" of the object households
-households <- editMap(kario_camp)
+
# Interactively draw a shape on your map
+# Your points are saved in the object "households" (which is a list)
+# More specifically in the slot "finished" of the object households
+households <- editMap(kario_map)
 
-## access your points 
+# Access your points 
 households <- households$finished

 
-

-
Task B: Read in existing rooftop points (optional)

+

+
Task B: Automatically identify rooftops (optional)

 
This is an alternative step that replaces marking roofs in the
 previous task

-
The full roofs identification of the refugee camp was in fact carried
-out semi-automatically, manually validated and saved in a shapefile. As
-with the polygon above we can use the sf package to read in
-the points to R.

+
It is possible to identify each rooftop in the camp automatically,
+using an OpenStreetMap feature that identifies buildings from satellite
+imagery. We can do this in R with the package osmdata.

+
First, we need to define a boundary box within which all the building
+data will be retrieved. We can get the boundary box for the Kario camp
+polygon as below:

+
# Get a boundary box for Kario camp from the shapefile:
+kario_bbox <- sf::st_bbox(kario_polygon)

+
Now we can retrieve the geometries (polygons and points) for all the
+buildings that fall within this boundary box. We use the
+opq function to query OpenStreetMap and select “building”
+as the feature that we want to retrieve:

+
# Get the builing points within the boundary box:
+buildings <- osmdata::opq(bbox = kario_bbox) %>%
+  
+  # Get the buildings in this boundary box:
+  osmdata::add_osm_feature(key = "building") %>%
+  
+  # Convert the building points to an osmdata_sf object:
+  osmdata::osmdata_sf() 

+
If you inspect the contents of buildings you can see
+that it is a list of different objects, including two data.frames:
+osm_points and osm_polygons. The coordinates
+in osm_points define the polygons for each building. Note
+that this means there are multiple points per building (one point for
+each corner in the shape that makes up the building). We can already see
+the outlines of the buildings on the map tiles however, so what we
+really need are a single pair of coordinates for each building (to
+represent one household). We can get this by calculating the centroid
+for each building polygon:

+
# Select the osm_polygons data.frame:
+households_osm <- buildings$osm_polygons %>% 
+  
+  # Convert it to a simple features object:
+  sf::st_as_sf() %>% 
+  
+  # Get centroids (points) for each building polygon:
+  mutate(points = sf::st_centroid(geometry))

+
Now we can add the points representing each household to the map:

+
# Update the map:
+kario_map <- kario_camp %>% 
+  
+  # Add the points representing buildings to the map:
+  addCircleMarkers(
+    
+    # Get building points to plot from the points column:
+    data = households_osm$points,  
+    
+    # Specify the radius of the circle markers:
+    radius = 5,                                          
+    
+    )
+
+
+# View the updated map:
+kario_map

+

+
+
Currently some of the buiding points are located outside the Kario
+camp polygon - this is because the query we used to fetch the building
+data used the boundary box of the polygon, not the polygon itself. We
+can remove the points outside the polygon by trimming the data with the
+st_intersection() function. This takes two arguments - x
+(the points data that we want to trim) and y (the polygon data that we
+want to use for the trimming).

+
households_osm <- households_osm %>% 
+
+  # Trim to points which are inside the Kario camp polygon:
+  sf::st_intersection(y = kario_polygon) %>% 
+  
+  # St_intersection adds an extra id column we don't need:
+  select(-id)

+
Now we can redraw the map and check it to make sure the extra points
+have been removed:

+
# Update the map of Kario camp:
+kario_map <- kario_camp %>% 
+  
+  # Plot the points representing buildings ontop:
+  addCircleMarkers(
+    data = households_osm$points,                    
+    radius = 5,             
+    )
+
+
+# View the map:
+kario_map

+

+
+
There is still a problem, however. The automated process has
+identified each building from satellite imagery, but not all of these
+buildings are occupied, and some of them are not households but have
+other functions (such as central storage facilities or bathroom blocks).
+In addition, some new families have arrived since the last satellite
+image was taken.

+
To ensure that each point represents an occupied household, the
+automatically generated points were validated by a field survey team and
+saved in a shapefile.

+
As with the polygon above we can use the sf package to
+read in the points to R.

 
An additional complication here is that the points are not stored as
 longitude and latitude, but as x and y coordinates. For this we use the
 st_transform() function to re-project the data to the
 appropriate format. We have to specify a coordinate reference system,
 and WGS84 is the standard one used by most mapping services online
 (e.g. Open Street Maps).

-
# read in your shapefile
-households <- read_sf(here::here("data", 
-                              "SDN_Kario_17Feb2019_roofs", 
-                              "SDN_Kario_17Feb2019_roofs.shp")) %>% 
-  # re-project x and y coordinates to be longitude and latitude. 
-  # note that crs 4326 is a code for WGS84
+
# Read in your shapefile
+households_validated <- sf::read_sf(
+  here::here("data", 
+             "SDN_Kario_17Feb2019_roofs",
+             "SDN_Kario_17Feb2019_roofs.shp")) %>% 
+  
+  # Re-project x and y coordinates to be longitude and latitude. 
+  # Note that crs 4326 is a code for WGS84
   st_transform(crs = 4326)

-
You could then view the polygon to your existing map.

-
kario_camp %>% 
-  # add circles for each point 
-  # (zoom in to see how they fit with rooftops)
-  addCircleMarkers(data = households)

-

-
+
Once again, we can redraw the map, this time using the manually
+validated points from the shapefile:

+
# Update the map of Kario camp:
+kario_map <- kario_camp %>% 
+  
+  # Plot the points representing buildings ontop:
+  addCircleMarkers(
+    data = households_validated,                    
+    radius = 5,             
+    )
+
+
+# View the map:
+kario_map

+

+
 

 

 
@@ -6781,7 +6892,7 @@ 
Task A: sampling points (T-Square)

 Kario survey, 400 roofs were sampled in order to do
 Roof-Simple-Random.

 
# select 50 random points from your polygon (i.e. not necessarily households)
-random_points <- st_sample(polygon, 50) %>% 
+random_points <- st_sample(kario_polygon, 50) %>% 
   # change back to a simple features dataframe (sampling made it a "collection")
   st_as_sf()
 
@@ -6790,9 +6901,9 @@ 
Task A: sampling points (T-Square)

 kario_camp %>% 
   # add circles for each point 
   # (zoom in to see how they fit with rooftops)
-  addCircleMarkers(data = random_points, color = "green")

-

-
+  addCircleMarkers(data = random_points, color = "green", radius = 5)
+
+

Task B: sampling roofs (Roof-Simple-Random)

@@ -6805,31 +6916,53 @@

Task B: sampling roofs (Roof-Simple-Random)

whereas sampling from points does not account for this. While it is not the case in this tutorial, in other situations you might have more than one polygon and need to sample roofs separately for each polygon.

-

Notice here when you visualise that the points are much closer to -where the population is actually situated in the camp (because you are -sampling households rather than any space within the camp boundary) - -this likely makes conducting your survey logistically simpler. In -addition, you only need to survey half as many houses compared to -T-Square (because you don’t need to find the next closest house).

-
# select 50 random points from your households
-random_roofs <- st_sample(households, 50) %>%
-  # transform back to points 
-  # (this comes out as a multipoint object, which is too complex)
+
We will first try to sample from the OSM household (unvalidated)
+data:

+
# Select 50 random roofs from the household_osm data:
+random_roofs_osm <- st_sample(households_osm, 50)
+
+# Visualise the selected points on a map:
+kario_camp %>% 
+  
+  # Add circles for each point 
+  addCircleMarkers(data = random_roofs_osm, 
+                   color = "orange", 
+                   radius = 5)

+

+
+
For comparison, we can also select 50 random points from the
+validated data. Note that this data requires some extra formatting
+before we can visualise it:

+
# Select 50 random points from your households
+random_roofs_val <- st_sample(households_validated, 50) %>%
+  
+  # Transform back to points 
+  # (this comes out as a multi-point object, which is too complex)
   st_cast("POINT") %>% 
-  # set the coordinate reference system for the points 
-  # note that crs 4326 is a code for WGS84
+  
+  # Set the coordinate reference system for the points 
+  # Note that crs 4326 is a code for WGS84
   st_set_crs(4326) %>% 
-  # change back to a simple features dataframe (multipoint made it a "collection")
+  
+  # Change back to a simple features data.frame (multi-point made it a "collection")
   st_as_sf()
   
 
 # Visualise the selected points 
 kario_camp %>% 
-  # add circles for each point 
-  # (zoom in to see how they fit with rooftops)
-  addCircleMarkers(data = random_roofs, color = "orange")

-

-
+  
+  # Add circles for each point 
+  addCircleMarkers(data = random_roofs_val, 
+                   color = "orange", 
+                   radius = 5)
+
+ +

Notice here when you visualise that the points are much closer to +where the population is actually situated in the camp (because you are +sampling households rather than any space within the camp boundary) - +this likely makes conducting your survey logistically simpler. In +addition, you only need to survey half as many houses compared to +T-Square (because you don’t need to find the next closest house).

@@ -6847,47 +6980,64 @@

Task A: Export the sampled points

only exported the random points (rather than the random roofs), but the same code can be used for either. To save our points we use the st_write() function from the sf package.

-
# export sampled points to gpx
+
# Export sampled points to gpx:
 random_points %>%
+  
   # gpx has a specific format which it accepts
-  # to meet this format we have to give a name to label each of the points
-  # note that this needs to be a lower case "name" 
-  # we consecutively number the points by combining "point" with the number of the row
+  # To meet this format we have to give a name to label each of the points
+  # Note that this needs to be a lower case "name" 
+  # We consecutively number the points by combining "point" with the number of the row
     mutate(name = paste0("point-", 1:n())) %>%
+  
+  # Now we can write the labelled points to a .gpx file:
     st_write(
-      # define the file path and extension (gpx)
-        dsn    = here::here("data", "random_points.gpx"),
-        # define which format to use (this could be guessed by the function)
+      
+      # Define the file path and extension (gpx)
+        dsn = here::here("data", "random_points.gpx"),
+        
+        # Define which format to use (this could be guessed by the function)
         driver = "GPX",
-        # overwrite any previously existing files. 
+        
+        # Overwrite any previously existing files. 
         delete_dsn = TRUE)

 
We could similarly use st_write() to save to a
 .shp or a .kml (used in Google Earth).

-
# export sampled points to shp 
+
# Export sampled points to shp 
 random_points %>%
-  # it can be useful to give a name to label each of the points
-  # we consecutively number the points by combining "point" with the number of the row
+  
+  # It can still be useful to give a name to label each of the points
+  # We consecutively number the points by combining "point" with the number of the row
     mutate(Name = paste0("point-", 1:n())) %>%
-    st_write(
-      # define the file path and extension (kml)
-        dsn    = here::here("data", "random_points.shp"),
-        # overwrite any previously existing files. 
+    
+  # Now we can write the labelled points to a shapefile:
+  st_write(
+      
+    # Define the file path and extension (kml)
+        dsn = here::here("data", "random_points.shp"),
+        
+        # Overwrite any previously existing files. 
         delete_dsn = TRUE)
 
 
-# export sampled points to kml 
+# Export sampled points to kml 
 random_points %>%
+  
   # kml has a specific format which it accepts
-  # to meet this format we have to give a name to label each of the points
-  # note that this needs to be a lower case "Name" 
-  # we consecutively number the points combining "point" with the number of the row
+  # To meet this format we have to give a name to label each of the points
+  # Note that this needs to be a lower case "Name" 
+  # We consecutively number the points combining "point" with the number of the row
     mutate(Name = paste0("point-", 1:n())) %>%
-    st_write(
-      # define the file path and extension (kml)
-        dsn    = here::here("data", "random_points.kml"),
-        # define which format to use (this could be guessed by the function)
+    
+  # Now we can write the labelled points to a kml file:
+  st_write(
+      
+    # Define the file path and extension (kml)
+        dsn = here::here("data", "random_points.kml"),
+        
+        # Define which format to use (this could be guessed by the function)
         driver = "kml",
-        # overwrite any previously existing files. 
+        
+        # Overwrite any previously existing files. 
         delete_dsn = TRUE)

 
 

@@ -6921,12 +7071,12 @@ 
Task B: Transfer the sampled points

 

 
 
-

-
Step 4: Calculate

+

+
Step 5: Calculate

 

 
Objective

 
The objective of this step is to calculate your population estimates
-using T-square and then Roof-Simple-Random.

+using T-square and then the Roof-Simple-Random method.

 

 
Task A: Enter T-Square data

 
We are now going to pretend that the field survey has been completed.
@@ -6989,7 +7139,7 @@ 
Task B: Calculate T-Square population estimate

 function form sf.

 
## calculate the area of within the polygon and return as a numeric
 ## otherwise it is returned as a measurement (with m2 on the end)
-camp_area <- st_area(polygon) %>% 
+camp_area <- st_area(kario_polygon) %>% 
   as.numeric()
 
 ## view your numbers
@@ -7027,7 +7177,7 @@ 
Task B: Calculate T-Square population estimate

 first_point <- survey_data %>% 
   select(gps_point) %>% 
   ## split in to a latitude and longitude column 
-  tidyr::separate(gps_point,into = c("lat","lon"), sep=", ") %>%
+  tidyr::separate(gps_point,into = c("lat","lon"), sep = ", ") %>%
   ## change columns to numeric
   mutate(across(c(lat, lon), as.numeric)) %>%
   ## set as sf points and define the crs
@@ -7036,14 +7186,14 @@ 
Task B: Calculate T-Square population estimate

 ## extract the first randomly sampled points
 second_point <- survey_data %>% 
   select(nearest_house_gps_point) %>% 
-  tidyr::separate(nearest_house_gps_point,into = c("lat","lon"), sep=", ") %>%
+  tidyr::separate(nearest_house_gps_point,into = c("lat","lon"), sep = ", ") %>%
   mutate(across(c(lat, lon), as.numeric)) %>%
   st_as_sf(coords = c("lat","lon"), crs = 4326)
 
 ## extract the first randomly sampled points
 third_point <- survey_data %>% 
   select(second_house_gps_point) %>% 
-  tidyr::separate(second_house_gps_point,into = c("lat","lon"), sep=", ") %>%
+  tidyr::separate(second_house_gps_point,into = c("lat","lon"), sep = ", ") %>%
   mutate(across(c(lat, lon), as.numeric)) %>%
   st_as_sf(coords = c("lat","lon"), crs = 4326)

 
Then we need to calculate the distances between our household GPS
@@ -7162,7 +7312,7 @@ 
Task D: Calculate Roof-Simple-Random population estimate

   summarise(across(c(under_five, five_plus, overall), mean))
 
 ## pull the number of households from the shapefile 
-num_hh <- nrow(households)
+num_hh <- nrow(households_validated)
 
 ## multiply the averages by the number of houses 
 population <- averages * num_hh
@@ -7197,7 +7347,7 @@ 
Appendix

 the boundary of an area using your mobile phone - OsmAnd then creates a
 .gpx file. You can read in these files with the following
 code.

-
## read in a track from osmand
+
## Read in a track from osmAnd
 polygon <- st_read(
   dsn = here::here("data", "track.gpx")
 )

diff --git a/docs/gis_sampling.html b/docs/gis_sampling.html
index 082e72d..7bab766 100644
--- a/docs/gis_sampling.html
+++ b/docs/gis_sampling.html
@@ -1667,20 +1667,18 @@
       document.body.style.height = "100%";
       document.documentElement.style.width = "100%";
       document.documentElement.style.height = "100%";
-      if (cel) {
-        cel.style.position = "absolute";
-        var pad = unpackPadding(sizing.padding);
-        cel.style.top = pad.top + "px";
-        cel.style.right = pad.right + "px";
-        cel.style.bottom = pad.bottom + "px";
-        cel.style.left = pad.left + "px";
-        el.style.width = "100%";
-        el.style.height = "100%";
-      }
+      cel.style.position = "absolute";
+      var pad = unpackPadding(sizing.padding);
+      cel.style.top = pad.top + "px";
+      cel.style.right = pad.right + "px";
+      cel.style.bottom = pad.bottom + "px";
+      cel.style.left = pad.left + "px";
+      el.style.width = "100%";
+      el.style.height = "100%";
 
       return {
-        getWidth: function() { return cel.offsetWidth; },
-        getHeight: function() { return cel.offsetHeight; }
+        getWidth: function() { return cel.getBoundingClientRect().width; },
+        getHeight: function() { return cel.getBoundingClientRect().height; }
       };
 
     } else {
@@ -1688,8 +1686,8 @@
       el.style.height = px(sizing.height);
 
       return {
-        getWidth: function() { return el.offsetWidth; },
-        getHeight: function() { return el.offsetHeight; }
+        getWidth: function() { return cel.getBoundingClientRect().width; },
+        getHeight: function() { return cel.getBoundingClientRect().height; }
       };
     }
   }
@@ -1913,8 +1911,8 @@
 
           elementData(el, "initialized", true);
           if (bindingDef.initialize) {
-            var result = bindingDef.initialize(el, el.offsetWidth,
-              el.offsetHeight);
+            var rect = el.getBoundingClientRect();
+            var result = bindingDef.initialize(el, rect.width, rect.height);
             elementData(el, "init_result", result);
           }
         }
@@ -1956,29 +1954,30 @@
       forEach(matches, function(el) {
         var sizeObj = initSizing(el, binding);
 
+        var getSize = function(el) {
+          if (sizeObj) {
+            return {w: sizeObj.getWidth(), h: sizeObj.getHeight()}
+          } else {
+            var rect = el.getBoundingClientRect();
+            return {w: rect.width, h: rect.height}
+          }
+        };
+
         if (hasClass(el, "html-widget-static-bound"))
           return;
         el.className = el.className + " html-widget-static-bound";
 
         var initResult;
         if (binding.initialize) {
-          initResult = binding.initialize(el,
-            sizeObj ? sizeObj.getWidth() : el.offsetWidth,
-            sizeObj ? sizeObj.getHeight() : el.offsetHeight
-          );
+          var size = getSize(el);
+          initResult = binding.initialize(el, size.w, size.h);
           elementData(el, "init_result", initResult);
         }
 
         if (binding.resize) {
-          var lastSize = {
-            w: sizeObj ? sizeObj.getWidth() : el.offsetWidth,
-            h: sizeObj ? sizeObj.getHeight() : el.offsetHeight
-          };
+          var lastSize = getSize(el);
           var resizeHandler = function(e) {
-            var size = {
-              w: sizeObj ? sizeObj.getWidth() : el.offsetWidth,
-              h: sizeObj ? sizeObj.getHeight() : el.offsetHeight
-            };
+            var size = getSize(el);
             if (size.w === 0 && size.h === 0)
               return;
             if (size.w === lastSize.w && size.h === lastSize.h)
@@ -2280,7 +2279,6 @@
     return result;
   }
 })();
-
 
 
 
+

+
 

 

 
Task B: Draw the polygon

@@ -6655,14 +6656,14 @@ 
Task B: Draw the polygon

 done, click on “Finish”. Then click on “Done” on the bottom right. Have
 a quick look at the mapedit
 website if you need additional help.

-
# interactively draw a shape on your map
-# your polygon is saved in the object "polygon" (which is a list)
-# more specifically in the slot "finished" of the object polygon.
-polygon <- editMap(kario_camp)
+
# Interactively draw a shape on your map
+# Your polygon is saved in the object "polygon" (which is a list)
+# More specifically in the slot "finished" of the object polygon.
+kario_polygon <- editMap(kario_camp)
 
 
-# access your polygon 
-polygon <- polygon$finished

+# Access your polygon 
+kario_polygon <- kario_polygon$finished

 

 

 
Task C: Read in an existing polygon (optional)

@@ -6678,24 +6679,28 @@ 
Task C: Read in an existing polygon (optional)

 a shapefile consists of many more files than just the .shp,
 if we tell R where that is, it then knows to pull in the others.

 
# read in your shapefile
-polygon <- read_sf(here::here("data", 
-                              "SDN_Kario_17Feb2019_polygon", 
-                              "SDN_Kario_17Feb2019_polygon.shp"))

+kario_polygon <- read_sf(here::here("data", 
+                                    "SDN_Kario_17Feb2019_polygon", 
+                                    "SDN_Kario_17Feb2019_polygon.shp"))

 
You could then view the polygon to your existing map.

 
kario_camp <- kario_camp %>% 
-  # plot your polygon on top of the existing tiles
+  
+  # Plot your polygon on top of the existing tiles
   addPolygons(
-    # specify your sf object
-        data  = polygon,
-        # color the border red
+    
+    # Specify your sf object
+        data  = kario_polygon,
+        
+        # Color the border red
         color = "red",
-        # make the middle see-through
+        
+        # Make the middle see-through
         fillOpacity = 0)
 
-# view your map
+# View the map:
 kario_camp

-

-
+

+
 
 
 
@@ -6723,42 +6728,148 @@ 
Task A: Mark the roofs

 “Finish”. Then click on “Done” on the bottom right. Have a quick look at
 the mapedit
 website if you need additional help.

-
# interactively draw a shape on your map
-# your points are saved in the object "households" (which is a list)
-# more specifically in the slot "finished" of the object households
-households <- editMap(kario_camp)
+
# Interactively draw a shape on your map
+# Your points are saved in the object "households" (which is a list)
+# More specifically in the slot "finished" of the object households
+households <- editMap(kario_map)
 
-## access your points 
+# Access your points 
 households <- households$finished

 
-

-
Task B: Read in existing rooftop points (optional)

+

+
Task B: Automatically identify rooftops (optional)

 
This is an alternative step that replaces marking roofs in the
 previous task

-
The full roofs identification of the refugee camp was in fact carried
-out semi-automatically, manually validated and saved in a shapefile. As
-with the polygon above we can use the sf package to read in
-the points to R.

+
It is possible to identify each rooftop in the camp automatically,
+using an OpenStreetMap feature that identifies buildings from satellite
+imagery. We can do this in R with the package osmdata.

+
First, we need to define a boundary box within which all the building
+data will be retrieved. We can get the boundary box for the Kario camp
+polygon as below:

+
# Get a boundary box for Kario camp from the shapefile:
+kario_bbox <- sf::st_bbox(kario_polygon)

+
Now we can retrieve the geometries (polygons and points) for all the
+buildings that fall within this boundary box. We use the
+opq function to query OpenStreetMap and select “building”
+as the feature that we want to retrieve:

+
# Get the builing points within the boundary box:
+buildings <- osmdata::opq(bbox = kario_bbox) %>%
+  
+  # Get the buildings in this boundary box:
+  osmdata::add_osm_feature(key = "building") %>%
+  
+  # Convert the building points to an osmdata_sf object:
+  osmdata::osmdata_sf() 

+
If you inspect the contents of buildings you can see
+that it is a list of different objects, including two data.frames:
+osm_points and osm_polygons. The coordinates
+in osm_points define the polygons for each building. Note
+that this means there are multiple points per building (one point for
+each corner in the shape that makes up the building). We can already see
+the outlines of the buildings on the map tiles however, so what we
+really need are a single pair of coordinates for each building (to
+represent one household). We can get this by calculating the centroid
+for each building polygon:

+
# Select the osm_polygons data.frame:
+households_osm <- buildings$osm_polygons %>% 
+  
+  # Convert it to a simple features object:
+  sf::st_as_sf() %>% 
+  
+  # Get centroids (points) for each building polygon:
+  mutate(points = sf::st_centroid(geometry))

+
Now we can add the points representing each household to the map:

+
# Update the map:
+kario_map <- kario_camp %>% 
+  
+  # Add the points representing buildings to the map:
+  addCircleMarkers(
+    
+    # Get building points to plot from the points column:
+    data = households_osm$points,  
+    
+    # Specify the radius of the circle markers:
+    radius = 5,                                          
+    
+    )
+
+
+# View the updated map:
+kario_map

+

+
+
Currently some of the buiding points are located outside the Kario
+camp polygon - this is because the query we used to fetch the building
+data used the boundary box of the polygon, not the polygon itself. We
+can remove the points outside the polygon by trimming the data with the
+st_intersection() function. This takes two arguments - x
+(the points data that we want to trim) and y (the polygon data that we
+want to use for the trimming).

+
households_osm <- households_osm %>% 
+
+  # Trim to points which are inside the Kario camp polygon:
+  sf::st_intersection(y = kario_polygon) %>% 
+  
+  # St_intersection adds an extra id column we don't need:
+  select(-id)

+
Now we can redraw the map and check it to make sure the extra points
+have been removed:

+
# Update the map of Kario camp:
+kario_map <- kario_camp %>% 
+  
+  # Plot the points representing buildings ontop:
+  addCircleMarkers(
+    data = households_osm$points,                    
+    radius = 5,             
+    )
+
+
+# View the map:
+kario_map

+

+
+
There is still a problem, however. The automated process has
+identified each building from satellite imagery, but not all of these
+buildings are occupied, and some of them are not households but have
+other functions (such as central storage facilities or bathroom blocks).
+In addition, some new families have arrived since the last satellite
+image was taken.

+
To ensure that each point represents an occupied household, the
+automatically generated points were validated by a field survey team and
+saved in a shapefile.

+
As with the polygon above we can use the sf package to
+read in the points to R.

 
An additional complication here is that the points are not stored as
 longitude and latitude, but as x and y coordinates. For this we use the
 st_transform() function to re-project the data to the
 appropriate format. We have to specify a coordinate reference system,
 and WGS84 is the standard one used by most mapping services online
 (e.g. Open Street Maps).

-
# read in your shapefile
-households <- read_sf(here::here("data", 
-                              "SDN_Kario_17Feb2019_roofs", 
-                              "SDN_Kario_17Feb2019_roofs.shp")) %>% 
-  # re-project x and y coordinates to be longitude and latitude. 
-  # note that crs 4326 is a code for WGS84
+
# Read in your shapefile
+households_validated <- sf::read_sf(
+  here::here("data", 
+             "SDN_Kario_17Feb2019_roofs",
+             "SDN_Kario_17Feb2019_roofs.shp")) %>% 
+  
+  # Re-project x and y coordinates to be longitude and latitude. 
+  # Note that crs 4326 is a code for WGS84
   st_transform(crs = 4326)

-
You could then view the polygon to your existing map.

-
kario_camp %>% 
-  # add circles for each point 
-  # (zoom in to see how they fit with rooftops)
-  addCircleMarkers(data = households)

-

-
+
Once again, we can redraw the map, this time using the manually
+validated points from the shapefile:

+
# Update the map of Kario camp:
+kario_map <- kario_camp %>% 
+  
+  # Plot the points representing buildings ontop:
+  addCircleMarkers(
+    data = households_validated,                    
+    radius = 5,             
+    )
+
+
+# View the map:
+kario_map

+

+
 

 

 
@@ -6781,7 +6892,7 @@ 
Task A: sampling points (T-Square)

 Kario survey, 400 roofs were sampled in order to do
 Roof-Simple-Random.

 
# select 50 random points from your polygon (i.e. not necessarily households)
-random_points <- st_sample(polygon, 50) %>% 
+random_points <- st_sample(kario_polygon, 50) %>% 
   # change back to a simple features dataframe (sampling made it a "collection")
   st_as_sf()
 
@@ -6790,9 +6901,9 @@ 
Task A: sampling points (T-Square)

 kario_camp %>% 
   # add circles for each point 
   # (zoom in to see how they fit with rooftops)
-  addCircleMarkers(data = random_points, color = "green")

-

-
+  addCircleMarkers(data = random_points, color = "green", radius = 5)

+

+
 
 

 
Task B: sampling roofs (Roof-Simple-Random)

@@ -6805,31 +6916,53 @@ 
Task B: sampling roofs (Roof-Simple-Random)

 whereas sampling from points does not account for this. While it is not
 the case in this tutorial, in other situations you might have more than
 one polygon and need to sample roofs separately for each polygon.

-
Notice here when you visualise that the points are much closer to
-where the population is actually situated in the camp (because you are
-sampling households rather than any space within the camp boundary) -
-this likely makes conducting your survey logistically simpler. In
-addition, you only need to survey half as many houses compared to
-T-Square (because you don’t need to find the next closest house).

-
# select 50 random points from your households
-random_roofs <- st_sample(households, 50) %>%
-  # transform back to points 
-  # (this comes out as a multipoint object, which is too complex)
+
We will first try to sample from the OSM household (unvalidated)
+data:

+
# Select 50 random roofs from the household_osm data:
+random_roofs_osm <- st_sample(households_osm, 50)
+
+# Visualise the selected points on a map:
+kario_camp %>% 
+  
+  # Add circles for each point 
+  addCircleMarkers(data = random_roofs_osm, 
+                   color = "orange", 
+                   radius = 5)

+

+
+
For comparison, we can also select 50 random points from the
+validated data. Note that this data requires some extra formatting
+before we can visualise it:

+
# Select 50 random points from your households
+random_roofs_val <- st_sample(households_validated, 50) %>%
+  
+  # Transform back to points 
+  # (this comes out as a multi-point object, which is too complex)
   st_cast("POINT") %>% 
-  # set the coordinate reference system for the points 
-  # note that crs 4326 is a code for WGS84
+  
+  # Set the coordinate reference system for the points 
+  # Note that crs 4326 is a code for WGS84
   st_set_crs(4326) %>% 
-  # change back to a simple features dataframe (multipoint made it a "collection")
+  
+  # Change back to a simple features data.frame (multi-point made it a "collection")
   st_as_sf()
   
 
 # Visualise the selected points 
 kario_camp %>% 
-  # add circles for each point 
-  # (zoom in to see how they fit with rooftops)
-  addCircleMarkers(data = random_roofs, color = "orange")

-

-
+  
+  # Add circles for each point 
+  addCircleMarkers(data = random_roofs_val, 
+                   color = "orange", 
+                   radius = 5)

+

+
+
Notice here when you visualise that the points are much closer to
+where the population is actually situated in the camp (because you are
+sampling households rather than any space within the camp boundary) -
+this likely makes conducting your survey logistically simpler. In
+addition, you only need to survey half as many houses compared to
+T-Square (because you don’t need to find the next closest house).

 

 
 
@@ -6847,47 +6980,64 @@ 
Task A: Export the sampled points

 only exported the random points (rather than the random roofs), but the
 same code can be used for either. To save our points we use the
 st_write() function from the sf package.

-
# export sampled points to gpx
+
# Export sampled points to gpx:
 random_points %>%
+  
   # gpx has a specific format which it accepts
-  # to meet this format we have to give a name to label each of the points
-  # note that this needs to be a lower case "name" 
-  # we consecutively number the points by combining "point" with the number of the row
+  # To meet this format we have to give a name to label each of the points
+  # Note that this needs to be a lower case "name" 
+  # We consecutively number the points by combining "point" with the number of the row
     mutate(name = paste0("point-", 1:n())) %>%
+  
+  # Now we can write the labelled points to a .gpx file:
     st_write(
-      # define the file path and extension (gpx)
-        dsn    = here::here("data", "random_points.gpx"),
-        # define which format to use (this could be guessed by the function)
+      
+      # Define the file path and extension (gpx)
+        dsn = here::here("data", "random_points.gpx"),
+        
+        # Define which format to use (this could be guessed by the function)
         driver = "GPX",
-        # overwrite any previously existing files. 
+        
+        # Overwrite any previously existing files. 
         delete_dsn = TRUE)

 
We could similarly use st_write() to save to a
 .shp or a .kml (used in Google Earth).

-
# export sampled points to shp 
+
# Export sampled points to shp 
 random_points %>%
-  # it can be useful to give a name to label each of the points
-  # we consecutively number the points by combining "point" with the number of the row
+  
+  # It can still be useful to give a name to label each of the points
+  # We consecutively number the points by combining "point" with the number of the row
     mutate(Name = paste0("point-", 1:n())) %>%
-    st_write(
-      # define the file path and extension (kml)
-        dsn    = here::here("data", "random_points.shp"),
-        # overwrite any previously existing files. 
+    
+  # Now we can write the labelled points to a shapefile:
+  st_write(
+      
+    # Define the file path and extension (kml)
+        dsn = here::here("data", "random_points.shp"),
+        
+        # Overwrite any previously existing files. 
         delete_dsn = TRUE)
 
 
-# export sampled points to kml 
+# Export sampled points to kml 
 random_points %>%
+  
   # kml has a specific format which it accepts
-  # to meet this format we have to give a name to label each of the points
-  # note that this needs to be a lower case "Name" 
-  # we consecutively number the points combining "point" with the number of the row
+  # To meet this format we have to give a name to label each of the points
+  # Note that this needs to be a lower case "Name" 
+  # We consecutively number the points combining "point" with the number of the row
     mutate(Name = paste0("point-", 1:n())) %>%
-    st_write(
-      # define the file path and extension (kml)
-        dsn    = here::here("data", "random_points.kml"),
-        # define which format to use (this could be guessed by the function)
+    
+  # Now we can write the labelled points to a kml file:
+  st_write(
+      
+    # Define the file path and extension (kml)
+        dsn = here::here("data", "random_points.kml"),
+        
+        # Define which format to use (this could be guessed by the function)
         driver = "kml",
-        # overwrite any previously existing files. 
+        
+        # Overwrite any previously existing files. 
         delete_dsn = TRUE)

 
 

@@ -6921,12 +7071,12 @@ 
Task B: Transfer the sampled points

 

 
 
-

-
Step 4: Calculate

+

+
Step 5: Calculate

 

 
Objective

 
The objective of this step is to calculate your population estimates
-using T-square and then Roof-Simple-Random.

+using T-square and then the Roof-Simple-Random method.

 

 
Task A: Enter T-Square data

 
We are now going to pretend that the field survey has been completed.
@@ -6989,7 +7139,7 @@ 
Task B: Calculate T-Square population estimate

 function form sf.

 
## calculate the area of within the polygon and return as a numeric
 ## otherwise it is returned as a measurement (with m2 on the end)
-camp_area <- st_area(polygon) %>% 
+camp_area <- st_area(kario_polygon) %>% 
   as.numeric()
 
 ## view your numbers
@@ -7027,7 +7177,7 @@ 
Task B: Calculate T-Square population estimate

 first_point <- survey_data %>% 
   select(gps_point) %>% 
   ## split in to a latitude and longitude column 
-  tidyr::separate(gps_point,into = c("lat","lon"), sep=", ") %>%
+  tidyr::separate(gps_point,into = c("lat","lon"), sep = ", ") %>%
   ## change columns to numeric
   mutate(across(c(lat, lon), as.numeric)) %>%
   ## set as sf points and define the crs
@@ -7036,14 +7186,14 @@ 
Task B: Calculate T-Square population estimate

 ## extract the first randomly sampled points
 second_point <- survey_data %>% 
   select(nearest_house_gps_point) %>% 
-  tidyr::separate(nearest_house_gps_point,into = c("lat","lon"), sep=", ") %>%
+  tidyr::separate(nearest_house_gps_point,into = c("lat","lon"), sep = ", ") %>%
   mutate(across(c(lat, lon), as.numeric)) %>%
   st_as_sf(coords = c("lat","lon"), crs = 4326)
 
 ## extract the first randomly sampled points
 third_point <- survey_data %>% 
   select(second_house_gps_point) %>% 
-  tidyr::separate(second_house_gps_point,into = c("lat","lon"), sep=", ") %>%
+  tidyr::separate(second_house_gps_point,into = c("lat","lon"), sep = ", ") %>%
   mutate(across(c(lat, lon), as.numeric)) %>%
   st_as_sf(coords = c("lat","lon"), crs = 4326)

 
Then we need to calculate the distances between our household GPS
@@ -7162,7 +7312,7 @@ 
Task D: Calculate Roof-Simple-Random population estimate

   summarise(across(c(under_five, five_plus, overall), mean))
 
 ## pull the number of households from the shapefile 
-num_hh <- nrow(households)
+num_hh <- nrow(households_validated)
 
 ## multiply the averages by the number of houses 
 population <- averages * num_hh
@@ -7197,7 +7347,7 @@ 
Appendix

 the boundary of an area using your mobile phone - OsmAnd then creates a
 .gpx file. You can read in these files with the following
 code.

-
## read in a track from osmand
+
## Read in a track from osmAnd
 polygon <- st_read(
   dsn = here::here("data", "track.gpx")
 )