Load Packages

task 1:

Read the NYC postal areas in Shapefiles into sf objects.

# read in a Shapefile. With the 'shp' file extension, sf knows it is a Shapefile.
nyc_sf <- st_read("R-spatial/Data/nyc/nyc_acs_tracts.shp")
## Reading layer `nyc_acs_tracts' from data source 
##   `/Users/khadijajallow/Documents/R/R-spatial/Data/nyc/nyc_acs_tracts.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 2166 features and 113 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.25559 ymin: 40.49612 xmax: -73.70001 ymax: 40.91553
## Geodetic CRS:  NAD83

task 2:

Read and process the NYS health facilities spreadsheet data. Create sf objects from geographic coordinates.

#tuning Health Facilities into spatial object

#read csv 
health_df <- read_csv("R-spatial/Data/NYS_Health_Facility.csv", 
                               show_col_types = FALSE, 
                               lazy = FALSE)
#rename latitude and longitude
health_df <- health_df %>%
  rename("lat"= `Facility Latitude`, "long" = `Facility Longitude`) 

#filter for specific lat/long values 
health_df <- health_df %>%
  filter(!is.na(lat),lat >= 40, lat <= 46,
         !is.na(long), long >= -80, long <= -72) 

#see df
#str(health_df)
#make into spatial object 
health_sf <- st_as_sf(health_df, 
                               coords = c("long", "lat"))
str(health_df)
## spc_tbl_ [3,843 × 36] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
##  $ Facility ID                 : num [1:3843] 204 620 1156 2589 3455 ...
##  $ Facility Name               : chr [1:3843] "Hospice at Lourdes" "Charles T Sitrin Health Care Center Inc" "East Side Nursing Home" "Wellsville Manor Care Center" ...
##  $ Short Description           : chr [1:3843] "HSPC" "NH" "NH" "NH" ...
##  $ Description                 : chr [1:3843] "Hospice" "Residential Health Care Facility - SNF" "Residential Health Care Facility - SNF" "Residential Health Care Facility - SNF" ...
##  $ Facility Open Date          : chr [1:3843] "06/01/1985" "02/01/1989" "08/01/1979" "02/01/1989" ...
##  $ Facility Address 1          : chr [1:3843] "4102 Old Vestal Road" "2050 Tilden Avenue" "62 Prospect St" "4192A Bolivar Road" ...
##  $ Facility Address 2          : chr [1:3843] NA NA NA NA ...
##  $ Facility City               : chr [1:3843] "Vestal" "New Hartford" "Warsaw" "Wellsville" ...
##  $ Facility State              : chr [1:3843] "New York" "New York" "New York" "New York" ...
##  $ Facility Zip Code           : chr [1:3843] "13850" "13413" "14569" "14895" ...
##  $ Facility Phone Number       : num [1:3843] 6.08e+09 3.16e+09 5.86e+09 5.86e+09 7.17e+09 ...
##  $ Facility Fax Number         : num [1:3843] NA NA NA NA NA ...
##  $ Facility Website            : chr [1:3843] NA NA NA NA ...
##  $ Facility County Code        : num [1:3843] 3 32 60 2 14 ...
##  $ Facility County             : chr [1:3843] "Broome" "Oneida" "Wyoming" "Allegany" ...
##  $ Regional Office ID          : num [1:3843] 3 3 1 1 1 7 1 7 5 7 ...
##  $ Regional Office             : chr [1:3843] "Central New York Regional Office" "Central New York Regional Office" "Western Regional Office - Buffalo" "Western Regional Office - Buffalo" ...
##  $ Main Site Name              : chr [1:3843] NA NA NA NA ...
##  $ Main Site Facility ID       : num [1:3843] NA NA NA NA NA ...
##  $ Operating Certificate Number: chr [1:3843] "0301501F" "3227304N" "6027303N" "0228305N" ...
##  $ Operator Name               : chr [1:3843] "Our Lady of Lourdes Memorial Hospital Inc" "Charles T Sitrin Health Care Center, Inc" "East Side Nursing Home Inc" "Wellsville Manor LLC" ...
##  $ Operator Address 1          : chr [1:3843] "169 Riverside Drive" "Box 1000 Tilden Avenue" "62 Prospect Street" "4192a Bolivar Road" ...
##  $ Operator Address 2          : chr [1:3843] NA NA NA NA ...
##  $ Operator City               : chr [1:3843] "Binghamton" "New Hartford" "Warsaw" "Wellsville" ...
##  $ Operator State              : chr [1:3843] "New York" "New York" "New York" "New York" ...
##  $ Operator Zip Code           : chr [1:3843] "13905" "13413" "14569" "14897" ...
##  $ Cooperator Name             : chr [1:3843] NA NA NA NA ...
##  $ Cooperator Address          : chr [1:3843] NA NA NA NA ...
##  $ Cooperator Address 2        : chr [1:3843] NA NA NA NA ...
##  $ Cooperator City             : chr [1:3843] NA NA NA NA ...
##  $ Cooperator State            : chr [1:3843] "New York" "New York" "New York" "New York" ...
##  $ Cooperator Zip Code         : chr [1:3843] NA NA NA NA ...
##  $ Ownership Type              : chr [1:3843] "Not for Profit Corporation" "Not for Profit Corporation" "Business Corporation" "LLC" ...
##  $ lat                         : num [1:3843] 42.1 43.1 42.7 42.1 43 ...
##  $ long                        : num [1:3843] -76 -75.2 -78.1 -78 -78.7 ...
##  $ Facility Location           : chr [1:3843] "(42.097095, -75.975243)" "(43.05497, -75.228828)" "(42.738979, -78.12867)" "(42.126461, -77.967834)" ...
##  - attr(*, "spec")=
##   .. cols(
##   ..   `Facility ID` = col_double(),
##   ..   `Facility Name` = col_character(),
##   ..   `Short Description` = col_character(),
##   ..   Description = col_character(),
##   ..   `Facility Open Date` = col_character(),
##   ..   `Facility Address 1` = col_character(),
##   ..   `Facility Address 2` = col_character(),
##   ..   `Facility City` = col_character(),
##   ..   `Facility State` = col_character(),
##   ..   `Facility Zip Code` = col_character(),
##   ..   `Facility Phone Number` = col_double(),
##   ..   `Facility Fax Number` = col_double(),
##   ..   `Facility Website` = col_character(),
##   ..   `Facility County Code` = col_double(),
##   ..   `Facility County` = col_character(),
##   ..   `Regional Office ID` = col_double(),
##   ..   `Regional Office` = col_character(),
##   ..   `Main Site Name` = col_character(),
##   ..   `Main Site Facility ID` = col_double(),
##   ..   `Operating Certificate Number` = col_character(),
##   ..   `Operator Name` = col_character(),
##   ..   `Operator Address 1` = col_character(),
##   ..   `Operator Address 2` = col_character(),
##   ..   `Operator City` = col_character(),
##   ..   `Operator State` = col_character(),
##   ..   `Operator Zip Code` = col_character(),
##   ..   `Cooperator Name` = col_character(),
##   ..   `Cooperator Address` = col_character(),
##   ..   `Cooperator Address 2` = col_character(),
##   ..   `Cooperator City` = col_character(),
##   ..   `Cooperator State` = col_character(),
##   ..   `Cooperator Zip Code` = col_character(),
##   ..   `Ownership Type` = col_character(),
##   ..   `Facility Latitude` = col_double(),
##   ..   `Facility Longitude` = col_double(),
##   ..   `Facility Location` = col_character()
##   .. )
##  - attr(*, "problems")=<externalptr>
st_crs(health_sf)
## Coordinate Reference System: NA
st_crs(health_sf) <- 4326 #  EPSG
st_crs(health_sf)
## Coordinate Reference System:
##   User input: EPSG:4326 
##   wkt:
## GEOGCRS["WGS 84",
##     ENSEMBLE["World Geodetic System 1984 ensemble",
##         MEMBER["World Geodetic System 1984 (Transit)"],
##         MEMBER["World Geodetic System 1984 (G730)"],
##         MEMBER["World Geodetic System 1984 (G873)"],
##         MEMBER["World Geodetic System 1984 (G1150)"],
##         MEMBER["World Geodetic System 1984 (G1674)"],
##         MEMBER["World Geodetic System 1984 (G1762)"],
##         MEMBER["World Geodetic System 1984 (G2139)"],
##         MEMBER["World Geodetic System 1984 (G2296)"],
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ENSEMBLEACCURACY[2.0]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["Horizontal component of 3D system."],
##         AREA["World."],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]
#filtering for specific zipcode and city to plot
data <-health_sf  %>% 
  dplyr::filter(!is.na(`Operator Zip Code`), !is.na(health_sf$`Facility City`), health_sf$`Facility City` == "Albany")

#plot 
ggplot(data) +
  geom_sf(aes(color=`Operator Zip Code`)) +
  coord_sf(xlim = c(-74.06, -73.85), default_crs = sf::st_crs(4326))

#same for brooklyn 
data <-health_sf  %>% 
  dplyr::filter(!is.na(`Operator Zip Code`), !is.na(health_sf$`Facility City`), health_sf$`Facility City` == "Brooklyn")

ggplot(data) +
  geom_sf(aes(color=`Facility County Code`)) +
  coord_sf(xlim = c(-74.06, -73.85), default_crs = sf::st_crs(4326))

Get NYC Retail Food Stores

#NYS Retail Food Stores 


retail_df <- read_csv("R-spatial/Data/nys_retail_food_store_xy.csv", 
                      show_col_types = FALSE, 
                      lazy = FALSE)
colnames(retail_df) <- c('county', 'Licnum', 'optype', 'etype', 'ename', 'dbname', 'snum', 'sname', 'alinetwo', 'alinethree', 'city', 'state', 'zip', 'sqft', 'location', 'coords', 'lat', 'long')


colnames(retail_df)
##  [1] "county"     "Licnum"     "optype"     "etype"      "ename"     
##  [6] "dbname"     "snum"       "sname"      "alinetwo"   "alinethree"
## [11] "city"       "state"      "zip"        "sqft"       "location"  
## [16] "coords"     "lat"        "long"
#select wanted values, remove na values, and filter for specific lat/long values 
retail_df<- retail_df %>%
  select(county, etype,dbname, city, zip,state,lat, long) %>%
  filter(!is.na(county), !is.na(etype), !is.na(city),!is.na(zip),!is.na(state), !is.na(lat),lat >= 40, lat <= 46,
         !is.na(long), long >= -80, long <= -72) 
#see df
#str(retail_df)

#make into spatial object 
retail_sf <- st_as_sf(retail_df, 
                      coords = c("long", "lat"))

st_crs(retail_sf)
## Coordinate Reference System: NA
st_crs(retail_sf) <- 4326 #  EPSG
st_crs(retail_sf)
## Coordinate Reference System:
##   User input: EPSG:4326 
##   wkt:
## GEOGCRS["WGS 84",
##     ENSEMBLE["World Geodetic System 1984 ensemble",
##         MEMBER["World Geodetic System 1984 (Transit)"],
##         MEMBER["World Geodetic System 1984 (G730)"],
##         MEMBER["World Geodetic System 1984 (G873)"],
##         MEMBER["World Geodetic System 1984 (G1150)"],
##         MEMBER["World Geodetic System 1984 (G1674)"],
##         MEMBER["World Geodetic System 1984 (G1762)"],
##         MEMBER["World Geodetic System 1984 (G2139)"],
##         MEMBER["World Geodetic System 1984 (G2296)"],
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ENSEMBLEACCURACY[2.0]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["Horizontal component of 3D system."],
##         AREA["World."],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]
#filtering for nyc counties 
city_data <-retail_sf  %>% 
  dplyr::filter(county %in% c("Bronx","Kings", "New York", "Richmond", "Queens"))

#plot 
ggplot(city_data) +
  geom_sf(aes(color= etype)) +
  coord_sf(xlim = c(-73.06, -74.7), default_crs = sf::st_crs(4326))

###Using Simple Mapping for Verification

mapview(city_data, zcol='etype', layer.name='Establishment Type')
mapview(data, zcol='Operator Zip Code', layer.name='Operator Zip Code')