Using sf package for spatial counting points in polygons via spatial join

Using sf package for spatial counting points in polygons…

Here are the steps for spatial counting points (i.e., Long and Lat coordinates) using the {sf} package

Note. To save the R code from this tutorial, please copy and paste to your RStudio the lines of code located in the gray boxes.

1. Set working directory (i.e.,Working Folder). In this folder, you need to place your data files and your R script file.

setwd("C:/MyRData/Using sf package for spatial counting points in polygons")

2. Upload libraries (i.e.,Packages). Always make sure to include the {tidiverse} package!
   Note. You will get warnings every time you run the packages for the first time.

library(sf)

## Linking to GEOS 3.9.1, GDAL 3.2.1, PROJ 7.2.1; sf_use_s2() is TRUE

library(tidyverse)

## -- Attaching packages --------------------------------------- tidyverse 1.3.2 --

## v ggplot2 3.4.2     v purrr   1.0.1
## v tibble  3.2.1     v dplyr   1.1.2
## v tidyr   1.3.0     v stringr 1.5.0
## v readr   2.1.3     v forcats 0.5.2
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

3. Upload the shape file from your working directory/folder.

polygons <- st_read("tractct_37800_0000_2010_s100_census_1_shp_wgs84.shp")

## Reading layer `tractct_37800_0000_2010_s100_census_1_shp_wgs84' from data source `C:\MyRData\Using sf package for spatial counting points in polygons\tractct_37800_0000_2010_s100_census_1_shp_wgs84.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 833 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -73.72778 ymin: 40.95095 xmax: -71.78724 ymax: 42.0506
## Geodetic CRS:  WGS 84

4. Look for coordinate reference systems (CRS) type in ‘Geodetic CRS:’.

polygons

## Simple feature collection with 833 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -73.72778 ymin: 40.95095 xmax: -71.78724 ymax: 42.0506
## Geodetic CRS:  WGS 84
## First 10 features:
##    STATEFP10 COUNTYFP10 TRACTCE10     GEOID10 NAME10        NAMELSAD10 MTFCC10
## 1         09        001    030300 09001030300    303  Census Tract 303   G5020
## 2         09        001    035300 09001035300    353  Census Tract 353   G5020
## 3         09        001    011000 09001011000    110  Census Tract 110   G5020
## 4         09        001    030100 09001030100    301  Census Tract 301   G5020
## 5         09        001    030200 09001030200    302  Census Tract 302   G5020
## 6         09        001    030400 09001030400    304  Census Tract 304   G5020
## 7         09        001    030500 09001030500    305  Census Tract 305   G5020
## 8         09        001    105100 09001105100   1051 Census Tract 1051   G5020
## 9         09        001    105200 09001105200   1052 Census Tract 1052   G5020
## 10        09        001    060100 09001060100    601  Census Tract 601   G5020
##    FUNCSTAT10  ALAND10 AWATER10  INTPTLAT10   INTPTLON10         GEOID_AFF1
## 1           S  8268946  5919424 +41.0563491 -073.4736510 14000US09001030300
## 2           S 13442741    46233 +41.1215292 -073.4784315 14000US09001035300
## 3           S  4449326  3865378 +41.0122609 -073.5766968 14000US09001011000
## 4           S 11024824    27898 +41.1006634 -073.4881986 14000US09001030100
## 5           S  4579450    45867 +41.0854982 -073.4644649 14000US09001030200
## 6           S  3552882        0 +41.0670034 -073.4899354 14000US09001030400
## 7           S  5351221     5422 +41.0808098 -073.5000669 14000US09001030500
## 8           S 31837542  1272596 +41.2338658 -073.3085559 14000US09001105100
## 9           S 39176156  1945200 +41.2806799 -073.2995681 14000US09001105200
## 10          S  3985964        0 +41.2140780 -073.2415497 14000US09001060100
##             GEOID_AFF2                       geometry
## 1  140000US09001030300 MULTIPOLYGON (((-73.46193 4...
## 2  140000US09001035300 MULTIPOLYGON (((-73.52206 4...
## 3  140000US09001011000 MULTIPOLYGON (((-73.57299 4...
## 4  140000US09001030100 MULTIPOLYGON (((-73.47458 4...
## 5  140000US09001030200 MULTIPOLYGON (((-73.47617 4...
## 6  140000US09001030400 MULTIPOLYGON (((-73.5091 41...
## 7  140000US09001030500 MULTIPOLYGON (((-73.49686 4...
## 8  140000US09001105100 MULTIPOLYGON (((-73.25958 4...
## 9  140000US09001105200 MULTIPOLYGON (((-73.24686 4...
## 10 140000US09001060100 MULTIPOLYGON (((-73.25108 4...

5. Read your Longitude and Latitude points from the .csv file.

points <- read_csv('Long_Lat_points.csv')

## Rows: 4294 Columns: 3
## -- Column specification --------------------------------------------------------
## Delimiter: ","
## chr (1): Incident Date
## dbl (2): Latitude, Longitude
## 
## i Use `spec()` to retrieve the full column specification for this data.
## i Specify the column types or set `show_col_types = FALSE` to quiet this message.

6. Convert the ‘Long_Lat_points.csv’ points to an sf object.

points_sf <- st_as_sf(points, coords = c("Longitude", "Latitude"))

7. Set CRS for the ‘points_sf’ file to be the same as the one from the ‘polygons’ shape file (i.e., WGS84).

st_crs(points_sf) <- st_crs(polygons)

8. Perform the spatial join between ‘points_sf’ and ‘polygons’ using st_join function.
   Note. Search ‘?st_join()’ for more information related to different ‘join=’ types.

points_sf_joined <- st_join(points_sf, polygons, join = st_within)

9. Create a new dataset from the ‘points_sf_joined’ file by counting the number of points within the census tracts.
   The newly created variable will be named ‘n’.

census_tract_points_count <- count(as_tibble(points_sf_joined), NAME10) %>%
                                   arrange(desc(n)) %>% print()

## # A tibble: 692 x 2
##    NAME10     n
##    <chr>  <int>
##  1 4171      99
##  2 3501      64
##  3 6905      44
##  4 5003      42
##  5 8006      40
##  6 4161      38
##  7 739       38
##  8 3508      37
##  9 5028      37
## 10 441       35
## # i 682 more rows

10. Perform the spatial join between ‘polygons’ file and ‘census_tract_points_count’ file using left_join function.

spatial_count_result <- left_join(polygons, census_tract_points_count, by = c("NAME10" = "NAME10"))

11. Quickly map the number of pizza deliveries by census tract using the {mapview} package.
    Display the number of pizza deliveries by census tract using the variable ‘n’.

library(mapview)
mapview(spatial_count_result, zcol = "n", layer.name = "Pizza Deliveries")

12. Finally, save the new ‘spatial_count_result’ shape file to the working directory.
    In your working directory, you will have four files named ‘spatial_count_result’ with file extensions: .dbf; .prj; .shp; and .shx

# Save the result as a shapefile (if needed)
st_write(spatial_count_result, "spatial_count_result.shp")

## Writing layer `spatial_count_result' to data source 
##   `spatial_count_result.shp' using driver `ESRI Shapefile'
## Writing 833 features with 15 fields and geometry type Multi Polygon.

Congratulations for completing all the 12 steps for creating a shape file with spatial counting points by census tract!

If you are new to R Programming Language, don’t give up. Your R skills will get better with time.