Week 7 Assignment

Loading packages:

library(tidyverse)
library(sf)
library(sp)
library(mapview)

Task 1:

lab_folder <- file.path(Sys.getenv("HOME"), "Downloads", "R-Spatial_I_Lab")
nyc_postal_areas <- st_read(file.path(lab_folder, "ZIP_CODE_040114 2", "ZIP_CODE_040114.shp")) 

## Reading layer `ZIP_CODE_040114' from data source 
##   `/Users/ayeshanaveed/Downloads/R-Spatial_I_Lab/ZIP_CODE_040114 2/ZIP_CODE_040114.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 263 features and 12 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 913129 ymin: 120020.9 xmax: 1067494 ymax: 272710.9
## Projected CRS: NAD83 / New York Long Island (ftUS)

Task 2:

health_facilities_file <- file.path(lab_folder, "NYS_Health_Facility.csv")
health_facilities <- read.csv(health_facilities_file)

# Remove the rows with missing coordinates
health_facilities <- health_facilities[!is.na(health_facilities$Facility.Longitude) & 
                                       !is.na(health_facilities$Facility.Latitude), ]

health_facilities_sf <- st_as_sf(health_facilities, coords = c("Facility.Longitude", "Facility.Latitude"), crs = 4326) 

Task 3:

library(readr)
retail_food_stores <- read_csv("/Users/ayeshanaveed/Downloads/R-Spatial_I_Lab/nys_retail_food_store_xy.csv")

## Rows: 29389 Columns: 18
## ── Column specification ──────────────────────────────────────────────────────
## Delimiter: ","
## chr (11): ï..County, Operation.Type, Establishment.Type, Entity.Name, DBA....
## dbl  (4): License.Number, Zip.Code, Y, X
## num  (1): Square.Footage
## lgl  (2): Address.Line.2, Address.Line.3
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

# Get the missing values in X and Y columns
sum(is.na(retail_food_stores$X))  

## [1] 5417

sum(is.na(retail_food_stores$Y)) 

## [1] 5417

# Filter out the rows with NA in X and Y
retail_food_stores <- retail_food_stores[!is.na(retail_food_stores$X) & !is.na(retail_food_stores$Y), ]


retail_food_stores_sf <- st_as_sf(retail_food_stores, coords = c("X", "Y"), crs = 4326)

Task 4:

library(mapview)
library(sf)

m1 <- mapview(nyc_postal_areas, 
              col.regions = "lightblue", 
              alpha.regions = 0.2, 
              layer.name = "NYC Postal Areas")

m2 <- mapview(health_facilities_sf, 
              col.regions = "pink", 
              cex = 3, 
              layer.name = "Health Facilities")

m3 <- mapview(retail_food_stores_sf, 
              col.regions = "green", 
              cex = 2, 
              layer.name = "Retail Food Stores")

m_combined <- m1 + m2 + m3
m_combined

Task 5:

st_write(nyc_postal_areas, "nyc_spatial_data.gpkg", layer = "postal_areas", driver = "GPKG", append = TRUE)

## Updating layer `postal_areas' to data source `nyc_spatial_data.gpkg' using driver `GPKG'
## Writing 263 features with 12 fields and geometry type Polygon.

st_write(health_facilities_sf, "nyc_spatial_data.gpkg", layer = "health_facilities", driver = "GPKG", append = TRUE)

## Updating layer `health_facilities' to data source `nyc_spatial_data.gpkg' using driver `GPKG'
## Writing 3848 features with 34 fields and geometry type Point.

st_write(retail_food_stores_sf, "nyc_spatial_data.gpkg", layer = "retail_food_stores", driver = "GPKG", append = TRUE)

## Updating layer `retail_food_stores' to data source `nyc_spatial_data.gpkg' using driver `GPKG'
## Writing 23972 features with 16 fields and geometry type Point.