This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.
Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.
#setwd("Your folder path with the Lab 2 data")
setwd("C:/Users/ssrini06/Box/Tufts/UEP236_SpatStat/LabExercises/Lab2")
getwd()
[1] "C:/Users/ssrini06/Box/Tufts/UEP236_SpatStat/LabExercises/Lab2"
#install libraries once only
#install.packages("sf")
#install.packages("terra")
#install.packages("exactextractr")
#install.packages("tmap")
#load libraries
library(terra)
library(sf)
library(tmap)
library(exactextractr)
SECTION 2.1 Reading shapefiles into terra and shapefiles into sf
# Read layer as a terra spatial object
nyc <- vect("nyc_neighborhood.shp")
#vector
nyc
class : SpatVector
geometry : polygons
dimensions : 195, 7 (geometries, attributes)
extent : -74.25559, -73.70001, 40.49612, 40.91553 (xmin, xmax, ymin, ymax)
source : nyc_neighborhood.shp
coord. ref. : lon/lat WGS84(DD)
#summary of the table
summary(nyc)
county_fip shape_area shape_leng
Length:195 Min. : 5573902 Min. : 11000
Class :character 1st Qu.: 19392084 1st Qu.: 23824
Mode :character Median : 32629789 Median : 30550
Mean : 43226938 Mean : 42012
3rd Qu.: 50237459 3rd Qu.: 41877
Max. :327756690 Max. :490427
ntacode boro_code ntaname
Length:195 Min. :1 Length:195
Class :character 1st Qu.:2 Class :character
Mode :character Median :3 Mode :character
Mean :3
3rd Qu.:4
Max. :5
boro_name
Length:195
Class :character
Mode :character
#attributes in the table
names(nyc)
[1] "county_fip" "shape_area" "shape_leng" "ntacode"
[5] "boro_code" "ntaname" "boro_name"
#table attribute names
head(as.data.frame(nyc))
#saving the attribute table in a dataframe
nyc.df <- as.data.frame(nyc)
#class
class(nyc)
[1] "SpatVector"
attr(,"package")
[1] "terra"
#map by area of polygon
plot(nyc, "shape_area")
#write the vector back to a shapefile
writeVector(nyc, "nyc_new.shp", overwrite=TRUE)
#read using sf to get a different type of spatial object
nyc_sf <- read_sf("nyc_neighborhood.shp")
#note that this is a different object than the one from terra
nyc_sf
Simple feature collection with 195 features and 7 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -74.25559 ymin: 40.49612 xmax: -73.70001 ymax: 40.91553
Geodetic CRS: WGS84(DD)
class(nyc_sf)
[1] "sf" "tbl_df" "tbl" "data.frame"
#same attributes
names(nyc_sf)
[1] "county_fip" "shape_area" "shape_leng" "ntacode"
[5] "boro_code" "ntaname" "boro_name" "geometry"
#plotting sf objects is slightly different
plot(st_geometry(nyc_sf))
plot(nyc_sf["shape_area"])
#writing to a shapefile
#write_sf(nyc_sf, "nyc_sf.shp", APPEND=F)
Section 2.2 Reading a raster
#Read raster
nyc_elev <- rast("be_NYC_025_agg30.tif")
#check the CRS for the raster and notice the epsg
nyc_elev
class : SpatRaster
dimensions : 834, 696, 1 (nrow, ncol, nlyr)
resolution : 30, 30 (x, y)
extent : 979137, 1000017, 2e+05, 225020 (xmin, xmax, ymin, ymax)
coord. ref. : NAD83 / New York Long Island (ftUS) (EPSG:2263)
source : be_NYC_025_agg30.tif
name : be_NYC_025_agg30
min value : -16.72868
max value : 141.94435
#class
class(nyc_elev)
[1] "SpatRaster"
attr(,"package")
[1] "terra"
#plot the raster
plot(nyc_elev)
#Saving a raster
writeRaster(nyc_elev, "nyc_elev.tif", overwrite=TRUE)
SECTION 2.3 Coordinate system, projection, etc (CRS)
# Coordinate system and projection
crs(nyc)
[1] "GEOGCRS[\"WGS84(DD)\",\n DATUM[\"WGS84\",\n ELLIPSOID[\"WGS84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1,\n ID[\"EPSG\",9001]]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"geodetic longitude\",east,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"geodetic latitude\",north,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]]]"
#this is unprojected using WGS84 Datum
#change the coordinate system to UTM zone 18 which is appropriate for NYC
newcrs_proj4string <- crs("+proj=utm +zone=18 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")
#Easier way of referring to the CRS using the epsg code
newcrs_epsg <- crs("+init=EPSG:32618")
newcrs_epsg<- crs("EPSG:32618")
#note that you don't have to include the package name terra::
# its helpful to know which package is used for a function
nyc_proj1 <- terra::project(nyc, newcrs_proj4string)
nyc_proj2 <- terra::project(nyc, newcrs_epsg)
#Both have the same projected coordinate system
nyc_proj1
class : SpatVector
geometry : polygons
dimensions : 195, 7 (geometries, attributes)
extent : 563069.7, 609762.3, 4483096, 4529952 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
nyc_proj2
class : SpatVector
geometry : polygons
dimensions : 195, 7 (geometries, attributes)
extent : 563069.7, 609762.3, 4483096, 4529952 (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618)
crs(nyc_proj1)
[1] "PROJCRS[\"unknown\",\n BASEGEOGCRS[\"unknown\",\n DATUM[\"World Geodetic System 1984\",\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]],\n ID[\"EPSG\",6326]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8901]]],\n CONVERSION[\"UTM zone 18N\",\n METHOD[\"Transverse Mercator\",\n ID[\"EPSG\",9807]],\n PARAMETER[\"Latitude of natural origin\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8801]],\n PARAMETER[\"Longitude of natural origin\",-75,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8802]],\n PARAMETER[\"Scale factor at natural origin\",0.9996,\n SCALEUNIT[\"unity\",1],\n ID[\"EPSG\",8805]],\n PARAMETER[\"False easting\",500000,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8806]],\n PARAMETER[\"False northing\",0,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8807]],\n ID[\"EPSG\",16018]],\n CS[Cartesian,2],\n AXIS[\"(E)\",east,\n ORDER[1],\n LENGTHUNIT[\"metre\",1,\n ID[\"EPSG\",9001]]],\n AXIS[\"(N)\",north,\n ORDER[2],\n LENGTHUNIT[\"metre\",1,\n ID[\"EPSG\",9001]]]]"
crs(nyc_proj2)
[1] "PROJCRS[\"WGS 84 / UTM zone 18N\",\n BASEGEOGCRS[\"WGS 84\",\n DATUM[\"World Geodetic System 1984\",\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4326]],\n CONVERSION[\"UTM zone 18N\",\n METHOD[\"Transverse Mercator\",\n ID[\"EPSG\",9807]],\n PARAMETER[\"Latitude of natural origin\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8801]],\n PARAMETER[\"Longitude of natural origin\",-75,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8802]],\n PARAMETER[\"Scale factor at natural origin\",0.9996,\n SCALEUNIT[\"unity\",1],\n ID[\"EPSG\",8805]],\n PARAMETER[\"False easting\",500000,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8806]],\n PARAMETER[\"False northing\",0,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8807]]],\n CS[Cartesian,2],\n AXIS[\"(E)\",east,\n ORDER[1],\n LENGTHUNIT[\"metre\",1]],\n AXIS[\"(N)\",north,\n ORDER[2],\n LENGTHUNIT[\"metre\",1]],\n USAGE[\n SCOPE[\"Engineering survey, topographic mapping.\"],\n AREA[\"Between 78°W and 72°W, northern hemisphere between equator and 84°N, onshore and offshore. Bahamas. Canada - Nunavut; Ontario; Quebec. Colombia. Cuba. Ecuador. Greenland. Haiti. Jamica. Panama. Turks and Caicos Islands. United States (USA). Venezuela.\"],\n BBOX[0,-78,84,-72]],\n ID[\"EPSG\",32618]]"
#notice that the coordinate system units change on the X and Y axes
plot(nyc_proj1)
# add=T can be used to add a layer to the existing plot
plot(nyc_proj2, "shape_area", add=T)
#notice that the projected layer will not plot because the two layers have different coordinate systems
plot(nyc)
plot(nyc_proj2, "shape_area", add=T)
#saving projected coordinate system information from a spatial object
NYC_proj_info <- crs(nyc_proj1)
NYC_proj_info
[1] "PROJCRS[\"unknown\",\n BASEGEOGCRS[\"unknown\",\n DATUM[\"World Geodetic System 1984\",\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]],\n ID[\"EPSG\",6326]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8901]]],\n CONVERSION[\"UTM zone 18N\",\n METHOD[\"Transverse Mercator\",\n ID[\"EPSG\",9807]],\n PARAMETER[\"Latitude of natural origin\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8801]],\n PARAMETER[\"Longitude of natural origin\",-75,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8802]],\n PARAMETER[\"Scale factor at natural origin\",0.9996,\n SCALEUNIT[\"unity\",1],\n ID[\"EPSG\",8805]],\n PARAMETER[\"False easting\",500000,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8806]],\n PARAMETER[\"False northing\",0,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8807]],\n ID[\"EPSG\",16018]],\n CS[Cartesian,2],\n AXIS[\"(E)\",east,\n ORDER[1],\n LENGTHUNIT[\"metre\",1,\n ID[\"EPSG\",9001]]],\n AXIS[\"(N)\",north,\n ORDER[2],\n LENGTHUNIT[\"metre\",1,\n ID[\"EPSG\",9001]]]]"
nyc_proj3 <- terra::project(nyc, NYC_proj_info)
nyc_proj3
class : SpatVector
geometry : polygons
dimensions : 195, 7 (geometries, attributes)
extent : 563069.7, 609762.3, 4483096, 4529952 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
plot(nyc_proj3, "shape_area")
#for sf objects there is a different projection tool
crs_epsg2263 <- st_crs(2263)
nyc_sf_proj <- sf::st_transform(nyc_sf, crs_epsg2263)
plot(nyc_elev)
plot(st_geometry(nyc_sf_proj), add=T)
#check CRS for a raster
st_crs(nyc_elev)
Coordinate Reference System:
User input: NAD83 / New York Long Island (ftUS)
wkt:
PROJCRS["NAD83 / New York Long Island (ftUS)",
BASEGEOGCRS["NAD83",
DATUM["North American Datum 1983",
ELLIPSOID["GRS 1980",6378137,298.257222101004,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4269]],
CONVERSION["Lambert Conic Conformal (2SP)",
METHOD["Lambert Conic Conformal (2SP)",
ID["EPSG",9802]],
PARAMETER["Latitude of false origin",40.1666666666667,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8821]],
PARAMETER["Longitude of false origin",-74,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8822]],
PARAMETER["Latitude of 1st standard parallel",41.0333333333333,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8823]],
PARAMETER["Latitude of 2nd standard parallel",40.6666666666667,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8824]],
PARAMETER["Easting at false origin",300000,
LENGTHUNIT["metre",1],
ID["EPSG",8826]],
PARAMETER["Northing at false origin",0,
LENGTHUNIT["metre",1],
ID["EPSG",8827]]],
CS[Cartesian,2],
AXIS["easting",east,
ORDER[1],
LENGTHUNIT["US survey foot",0.304800609601219]],
AXIS["northing",north,
ORDER[2],
LENGTHUNIT["US survey foot",0.304800609601219]],
ID["EPSG",2263]]
crs(nyc_elev)
[1] "PROJCRS[\"NAD83 / New York Long Island (ftUS)\",\n BASEGEOGCRS[\"NAD83\",\n DATUM[\"North American Datum 1983\",\n ELLIPSOID[\"GRS 1980\",6378137,298.257222101004,\n LENGTHUNIT[\"metre\",1]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n ID[\"EPSG\",4269]],\n CONVERSION[\"Lambert Conic Conformal (2SP)\",\n METHOD[\"Lambert Conic Conformal (2SP)\",\n ID[\"EPSG\",9802]],\n PARAMETER[\"Latitude of false origin\",40.1666666666667,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8821]],\n PARAMETER[\"Longitude of false origin\",-74,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8822]],\n PARAMETER[\"Latitude of 1st standard parallel\",41.0333333333333,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8823]],\n PARAMETER[\"Latitude of 2nd standard parallel\",40.6666666666667,\n ANGLEUNIT[\"degree\",0.0174532925199433],\n ID[\"EPSG\",8824]],\n PARAMETER[\"Easting at false origin\",300000,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8826]],\n PARAMETER[\"Northing at false origin\",0,\n LENGTHUNIT[\"metre\",1],\n ID[\"EPSG\",8827]]],\n CS[Cartesian,2],\n AXIS[\"easting\",east,\n ORDER[1],\n LENGTHUNIT[\"US survey foot\",0.304800609601219]],\n AXIS[\"northing\",north,\n ORDER[2],\n LENGTHUNIT[\"US survey foot\",0.304800609601219]],\n ID[\"EPSG\",2263]]"
#different ways to store the CRS of the raster
raster_st_crs <- crs(nyc_elev)
raster_crs <- crs("+init=EPSG:2263")
raster_crs <- crs("EPSG:2263")
nyc_proj3 <- terra::project(nyc, raster_crs)
SECTION 2.4 Geocoding Lat and Long coordinates
#plot tweets and notice this is just a graph not a spatial object yet
nyc_tweets <- read.csv(file="NYC_Tweets.csv")
plot(nyc_tweets$Lon, nyc_tweets$Lat, pch=16, cex=0.5, col="blue")
#convert to a spatial object
nyc_tweet_pts <- st_as_sf(nyc_tweets, coords = c("Lon", "Lat"), crs = 4326)
#notice that this a different class than nyc
#terra creates vectors that are not the same as sf
class(nyc_tweet_pts)
[1] "sf" "data.frame"
plot(nyc)
plot(nyc_tweet_pts, col="blue", add=T)
Warning: ignoring all but the first attribute
# you can convert it to a spatial vector using
nyc_tweet_pts_vect <- vect(nyc_tweet_pts)
class(nyc_tweet_pts_vect)
[1] "SpatVector"
attr(,"package")
[1] "terra"
#optionally you can convert spat vectors to sf
nyc_sf <- sf::st_as_sf(nyc)
class(nyc_sf)
[1] "sf" "data.frame"
#If you leave it as a sf object
NYC_proj_info <- crs(nyc_proj1)
nyc_tweet_pts_pr <- st_transform(nyc_tweet_pts, NYC_proj_info)
plot(nyc_proj1)
plot(nyc_tweet_pts_pr, add=T)
Warning: ignoring all but the first attribute
#write an sf object to a shapefile
#st_write(nyc_tweet_pts, "nyc_tweets.shp", append=FALSE)
# insert your code here to project the tweets to match the raster CRS
# plot the tweets so you can see them with the elevation
SECTION 3.1 Attribute Joins (Table joins with a key)
# attribute joins using keys
nyc_popntable <- read.csv(file="nyc_population_neighborhood.csv")
#ntacode is the common key in both tables
nyc_neighborhood_pop <- merge(nyc, nyc_popntable, by.x = "ntacode", by.y = "ntacode")
#note the class is a spatVector
class(nyc_neighborhood_pop)
[1] "SpatVector"
attr(,"package")
[1] "terra"
names(nyc_neighborhood_pop)
[1] "ntacode" "county_fip" "shape_area" "shape_leng"
[5] "boro_code" "ntaname" "boro_name" "Borough"
[9] "FIPSCounty" "NTA.Name" "Pop2000" "Pop2010"
plot(nyc_neighborhood_pop, "Pop2010")
#creating a new attribute for density
nyc_neighborhood_pop$density10 = (nyc_neighborhood_pop$Pop2010/nyc_neighborhood_pop$shape_area)*1000
plot(nyc_neighborhood_pop, "density10")
#needs to convert to sf to make maps in tmap
nyc_neighborhood_pop_sf <- st_as_sf(nyc_neighborhood_pop)
qtm(nyc_neighborhood_pop_sf, fill = "density10")
NA
NA
SECTION 3.2 Spatial join (overlaying two vectors)
#check the CRS of the objects you want to join
st_crs(nyc_neighborhood_pop)
Coordinate Reference System:
User input: WGS84(DD)
wkt:
GEOGCRS["WGS84(DD)",
DATUM["WGS84",
ELLIPSOID["WGS84",6378137,298.257223563,
LENGTHUNIT["metre",1,
ID["EPSG",9001]]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic longitude",east,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic latitude",north,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]]]
st_crs(nyc_tweet_pts_vect)
Coordinate Reference System:
User input: WGS 84
wkt:
GEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
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]]
#also make sure they are the same spatial class
class(nyc_tweet_pts)
[1] "sf" "data.frame"
class(nyc_tweet_pts_vect)
[1] "SpatVector"
attr(,"package")
[1] "terra"
class(nyc_neighborhood_pop)
[1] "SpatVector"
attr(,"package")
[1] "terra"
class(nyc_sf)
[1] "sf" "data.frame"
#Spatial Joins st_join(target_sf, source_sf)
#works with sf objects so convert them to sf if they are spatVector
tweets_sf <- sf::st_as_sf(nyc_tweet_pts)
nyc_popn_sf <- sf::st_as_sf(nyc_neighborhood_pop)
class(tweets_sf)
[1] "sf" "data.frame"
class(nyc_popn_sf)
[1] "sf" "data.frame"
# use the project tool to make them both projected to EPSG 2263
crs_epsg2263 <- crs("+init=EPSG:2263")
nyc_popn_epsg2263 <- st_transform(nyc_popn_sf, crs_epsg2263)
tweets_sf_2263 <- st_transform(tweets_sf, crs_epsg2263)
# Then run spatial join
nhood_target <- st_join(nyc_popn_epsg2263, tweets_sf_2263)
tweets_target <- st_join(tweets_sf_2263, nyc_popn_epsg2263)
# Notice that the first layer is the target layer
#you will get polygons if the target is polygons
# and points if the target is points
nhood_target
Simple feature collection with 12130 features and 19 fields
Geometry type: GEOMETRY
Dimension: XY
Bounding box: xmin: 913175.1 ymin: 120121.9 xmax: 1067383 ymax: 272844.3
Projected CRS: NAD83 / New York Long Island (ftUS)
First 10 features:
ntacode county_fip shape_area shape_leng boro_code
1 BK88 047 54005019 39247.23 3
1.1 BK88 047 54005019 39247.23 3
1.2 BK88 047 54005019 39247.23 3
1.3 BK88 047 54005019 39247.23 3
1.4 BK88 047 54005019 39247.23 3
1.5 BK88 047 54005019 39247.23 3
1.6 BK88 047 54005019 39247.23 3
1.7 BK88 047 54005019 39247.23 3
1.8 BK88 047 54005019 39247.23 3
1.9 BK88 047 54005019 39247.23 3
ntaname boro_name Borough FIPSCounty NTA.Name
1 Borough Park Brooklyn Brooklyn 47 Borough Park
1.1 Borough Park Brooklyn Brooklyn 47 Borough Park
1.2 Borough Park Brooklyn Brooklyn 47 Borough Park
1.3 Borough Park Brooklyn Brooklyn 47 Borough Park
1.4 Borough Park Brooklyn Brooklyn 47 Borough Park
1.5 Borough Park Brooklyn Brooklyn 47 Borough Park
1.6 Borough Park Brooklyn Brooklyn 47 Borough Park
1.7 Borough Park Brooklyn Brooklyn 47 Borough Park
1.8 Borough Park Brooklyn Brooklyn 47 Borough Park
1.9 Borough Park Brooklyn Brooklyn 47 Borough Park
Pop2000 Pop2010 density10 FID F1 Tweet_ID User_ID
1 101055 106357 1.969391 67 NA 6.834556e+17 163232855
1.1 101055 106357 1.969391 270 NA 6.834592e+17 15164780
1.2 101055 106357 1.969391 941 NA 6.834717e+17 229176278
1.3 101055 106357 1.969391 1047 NA 6.834733e+17 509024947
1.4 101055 106357 1.969391 1219 NA 6.834769e+17 138578431
1.5 101055 106357 1.969391 1222 NA 6.834769e+17 509024947
1.6 101055 106357 1.969391 1785 NA 6.834892e+17 15164780
1.7 101055 106357 1.969391 1879 NA 6.834915e+17 509024947
1.8 101055 106357 1.969391 4977 NA 6.836577e+17 30497066
1.9 101055 106357 1.969391 5177 NA 6.836636e+17 509024947
Content
1 #coneyisland @ Coney Island Beach https://t.co/6wGdkIxzQi
1.1 Beis Yaakov play on a Saturday night. True story. (@ Franklin D. Roosevelt High School in Brooklyn NY) https://t.co/V3FPbnhNZa
1.2 @beastieweenie emerging from the frigid water at Coney Island yesterday during the New Years Day… https://t.co/2rnGP4tfo6
1.3 CHECK OUT THE FULL VIDEO ON OUR YOUTUBE FACEBOOK AND VIMEO CHANNELS! Dark brown #ouroboros… https://t.co/yOKgU9K9GO
1.4 Judges of reality shows: when a transgender or an ugly weird contestant shows up quit saying omg he's so cute! When he's obviously not
1.5 Sweet #seaturtle #momanddaughter #tattoo by sandydex_tattoos @tat2wonderland #tattoowonderland… https://t.co/FVJSOfm64t
1.6 I'm at @Sprinkles in Brooklyn NY https://t.co/S3MBJYkulh
1.7 CHECK OUT THE FULL VIDEO ON OUR YOUTUBE FACEBOOK AND VIMEO CHANNELS! Sweet #seaturtle… https://t.co/gdQoGTSjGl
1.8 Metaphor #deep #2016 @ Ft Hamilton Brooklyn Ny https://t.co/pjvTh6rz1i
1.9 @msfin_ conjured a #darkmark #tattoo for her brother @tat2wonderland #tattoowonderland #brooklyn… https://t.co/sSn6KZvDZu
Timestamp geometry
1 2016-01-03 01:11:19 POLYGON ((990897.9 169268.1...
1.1 2016-01-03 01:25:45 POLYGON ((990897.9 169268.1...
1.2 2016-01-03 02:15:35 POLYGON ((990897.9 169268.1...
1.3 2016-01-03 02:21:49 POLYGON ((990897.9 169268.1...
1.4 2016-01-03 02:36:02 POLYGON ((990897.9 169268.1...
1.5 2016-01-03 02:36:05 POLYGON ((990897.9 169268.1...
1.6 2016-01-03 03:25:01 POLYGON ((990897.9 169268.1...
1.7 2016-01-03 03:34:04 POLYGON ((990897.9 169268.1...
1.8 2016-01-03 14:34:29 POLYGON ((990897.9 169268.1...
1.9 2016-01-03 14:58:05 POLYGON ((990897.9 169268.1...
tweets_target
Simple feature collection with 12125 features and 19 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 915502.5 ymin: 124780.4 xmax: 1065026 ymax: 269580.3
Projected CRS: NAD83 / New York Long Island (ftUS)
First 10 features:
FID F1 Tweet_ID User_ID
1 0 NA 6.834544e+17 73683589
2 1 NA 6.834544e+17 59572535
3 2 NA 6.834544e+17 112073675
4 3 NA 6.834544e+17 406617881
5 4 NA 6.834544e+17 1475591984
6 5 NA 6.834545e+17 68322711
7 6 NA 6.834545e+17 22204765
8 7 NA 6.834545e+17 93502745
9 8 NA 6.834545e+17 23258204
10 9 NA 6.834545e+17 36410884
Content
1 by @chiptography at I/O - Chip Music NYC new year 2016 @ Bushwhick Brooklyn https://t.co/JKJVCEgkJB
2 BEST TIME EVER.. Later 2015🖕ðŸ\u008f½âœŒðŸ\u008f½ï¸\u008f.. And 2016 WAS GUD?!? 😜ðŸ\u008d¾ðŸŽ‰âœ”ï¸\u008f💯â\u009d—ï¸\u008f Now… https://t.co/EgrSna7zxm
3 These lights will inspire you. #empirestatebuilding #nyc @ Gantry Plaza State Park https://t.co/1dO0Ebmzrr
4 From #light to #dark #art #culture #newyorkcity #newyork #nyc @ Whitney Museum of American Art https://t.co/Uw6nEE1xLa
5 a fond farewell to this kitty always giving shade #✌ðŸ\u008f½ï¸\u008f @ Hellcat HeadQuarters https://t.co/dPN0xy7r27
6 Back at the sewing machine it's about to go down!! #creative #sewingproject #makingart… https://t.co/AZ7eLNhppc
7 Pretty cool poster. Who this character is? #fukuplus #mapache #nyc #poster @ Má Pêche (Momofuku) https://t.co/XZWKrC24RP
8 I'm at Gulf Gas Station in Brooklyn NY https://t.co/aS4GCcqM6H
9 About to order the menu @ STK https://t.co/XGI2eGjSXp
10 tonite! #doubleheadeddisco presents POST. a #postpunk #postholiday hangover. @nowherenyc 10pm no… https://t.co/lRCBGGlHHZ
Timestamp ntacode county_fip shape_area
1 2016-01-03 01:06:40 BK77 047 24927927
2 2016-01-03 01:06:41 MN04 061 16093788
3 2016-01-03 01:06:41 QN31 081 102350779
4 2016-01-03 01:06:48 MN23 061 25000527
5 2016-01-03 01:06:48 BK95 047 14522604
6 2016-01-03 01:06:55 QN15 081 54160919
7 2016-01-03 01:07:01 MN17 061 30192057
8 2016-01-03 01:07:06 BK82 047 117083807
9 2016-01-03 01:07:14 MN17 061 30192057
10 2016-01-03 01:07:14 MN22 061 10896915
shape_leng boro_code ntaname
1 26321.63 3 Bushwick North
2 17410.82 1 Hamilton Heights
3 74605.80 4 Hunters Point-Sunnyside-West Maspeth
4 29385.03 1 West Village
5 18756.70 3 Erasmus
6 48676.73 4 Far Rockaway-Bayswater
7 27035.74 1 Midtown-Midtown South
8 89197.32 3 East New York
9 27035.74 1 Midtown-Midtown South
10 13539.25 1 East Village
boro_name Borough FIPSCounty
1 Brooklyn Brooklyn 47
2 Manhattan Manhattan 61
3 Queens Queens 81
4 Manhattan Manhattan 61
5 Brooklyn Brooklyn 47
6 Queens Queens 81
7 Manhattan Manhattan 61
8 Brooklyn Brooklyn 47
9 Manhattan Manhattan 61
10 Manhattan Manhattan 61
NTA.Name Pop2000 Pop2010
1 Bushwick North 56093 57138
2 Hamilton Heights 50555 48520
3 Hunters Point-Sunnyside-West Maspeth 61947 63271
4 West Village 68483 66880
5 Erasmus 31392 29938
6 Far Rockaway-Bayswater 48344 50058
7 Midtown-Midtown South 25807 28630
8 East New York 83275 91958
9 Midtown-Midtown South 25807 28630
10 East Village 41746 44136
density10 geometry
1 2.2921280 POINT (1006204 195083.9)
2 3.0148279 POINT (997658.3 239666.1)
3 0.6181780 POINT (995641.6 210804.1)
4 2.6751437 POINT (981951.8 208698.6)
5 2.0614761 POINT (998124.4 176092.9)
6 0.9242458 POINT (1053445 161107.9)
7 0.9482627 POINT (990977.9 217204.1)
8 0.7854032 POINT (1016411 181983.7)
9 0.9482627 POINT (988975.1 214452.4)
10 4.0503208 POINT (988349.4 205982.8)
#Aggregate by ntaname to map by polygon
#this could take some time so wait, meditate, sing a song
tweets_aggregated_nhood <- aggregate(x = nhood_target, by = list(nhood_target$ntaname), FUN = length)
#remove additional columns
tweets_aggregated_nhood <- tweets_aggregated_nhood[,1:2]
#rename columns
colnames(tweets_aggregated_nhood) <- c("ntaname", "count","geometry")
#map
plot(tweets_aggregated_nhood["count"])
#map using qtm looks nicer
qtm(tweets_aggregated_nhood, fill = "count")
#join the population table notice that the keys have different names in each table
nyc_neighborhood_pop_tweets <- merge(tweets_aggregated_nhood, nyc_popntable, by.x = "ntaname", by.y = "NTA.Name")
#calculate tweets per capita
nyc_neighborhood_pop_tweets$tweet_per_capita <- (nyc_neighborhood_pop_tweets$count/nyc_neighborhood_pop_tweets$Pop2010)*10000
qtm(nyc_neighborhood_pop_tweets, fill = "tweet_per_capita")
SECTION 3.3 Zonal statistics (Overalying a raster and a vector)
#Zonal stats
plot(nyc_elev)
plot(st_geometry(nyc_popn_epsg2263), add = T)
#extract works like zonal statistics
extract_elev_bynhood <- extract(nyc_elev, nyc_popn_epsg2263, na.rm=TRUE, fun=mean)
#the output is a table
plot(extract_elev_bynhood)
#faster tool for zonal statistics from the exactextractr library
#this lets us also get a key ntacode to do a join
#whoever came up with exactextractr? a typo waiting to happen
extract_elev_bynhood2 <- exact_extract(
x = nyc_elev, # raster
y = nyc_popn_epsg2263, # vector zones
fun = "mean",
append_cols = "ntacode")
|
| | 0%
|
| | 1%
|
|= | 1%
|
|= | 2%
|
|= | 3%
|
|== | 3%
|
|== | 4%
|
|== | 5%
|
|=== | 5%
|
|=== | 6%
|
|=== | 7%
|
|==== | 7%
|
|==== | 8%
|
|==== | 9%
|
|===== | 9%
|
|===== | 10%
|
|===== | 11%
|
|====== | 11%
|
|====== | 12%
|
|======= | 13%
|
|======= | 14%
|
|======== | 15%
|
|======== | 16%
|
|========= | 17%
|
|========= | 18%
|
|========== | 19%
|
|========== | 20%
|
|========== | 21%
|
|=========== | 21%
|
|=========== | 22%
|
|============ | 23%
|
|============ | 24%
|
|============= | 25%
|
|============= | 26%
|
|============== | 27%
|
|============== | 28%
|
|=============== | 29%
|
|=============== | 30%
|
|================ | 31%
|
|================ | 32%
|
|================= | 33%
|
|================= | 34%
|
|================== | 34%
|
|================== | 35%
|
|================== | 36%
|
|=================== | 36%
|
|=================== | 37%
|
|=================== | 38%
|
|==================== | 38%
|
|==================== | 39%
|
|==================== | 40%
|
|===================== | 41%
|
|===================== | 42%
|
|====================== | 43%
|
|====================== | 44%
|
|======================= | 45%
|
|======================= | 46%
|
|======================== | 46%
|
|======================== | 47%
|
|======================== | 48%
|
|========================= | 48%
|
|========================= | 49%
|
|========================= | 50%
|
|========================== | 50%
|
|========================== | 51%
|
|========================== | 52%
|
|=========================== | 52%
|
|=========================== | 53%
|
|=========================== | 54%
|
|============================ | 54%
|
|============================ | 55%
|
|============================= | 56%
|
|============================= | 57%
|
|============================== | 58%
|
|============================== | 59%
|
|=============================== | 60%
|
|=============================== | 61%
|
|=============================== | 62%
|
|================================ | 62%
|
|================================ | 63%
|
|================================ | 64%
|
|================================= | 64%
|
|================================= | 65%
|
|================================= | 66%
|
|================================== | 66%
|
|================================== | 67%
|
|=================================== | 68%
|
|=================================== | 69%
|
|==================================== | 70%
|
|==================================== | 71%
|
|===================================== | 72%
|
|===================================== | 73%
|
|====================================== | 74%
|
|====================================== | 75%
|
|======================================= | 76%
|
|======================================= | 77%
|
|======================================== | 78%
|
|======================================== | 79%
|
|========================================= | 79%
|
|========================================= | 80%
|
|========================================= | 81%
|
|========================================== | 82%
|
|========================================== | 83%
|
|=========================================== | 84%
|
|=========================================== | 85%
|
|============================================ | 86%
|
|============================================ | 87%
|
|============================================= | 88%
|
|============================================= | 89%
|
|============================================== | 89%
|
|============================================== | 90%
|
|============================================== | 91%
|
|=============================================== | 91%
|
|=============================================== | 92%
|
|=============================================== | 93%
|
|================================================ | 93%
|
|================================================ | 94%
|
|================================================ | 95%
|
|================================================= | 95%
|
|================================================= | 96%
|
|================================================= | 97%
|
|================================================== | 97%
|
|================================================== | 98%
|
|================================================== | 99%
|
|===================================================| 99%
|
|===================================================| 100%
#attribute join it back to the neighborhood boundaries using ntacode
nyc_neighborhood_avgelev <- merge(nyc_popn_epsg2263, extract_elev_bynhood2, by.x = "ntacode", by.y = "ntacode")
plot(nyc_neighborhood_avgelev["mean"])
qtm(nyc_neighborhood_avgelev, fill = "mean")