-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathcensus.R
61 lines (39 loc) · 1.88 KB
/
census.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
#Currently working on mapping functions in the is R code, with the intention of pulling in Census data.
library(tidycensus)
library(rgdal)
library(sf)
library(raster)
#TOKEN (CENSUS HERE)
#ckey="yourtoken"
census_api_key(ckey, install = TRUE)
url<-"https://www2.census.gov/programs-surveys/popest/geographies/2016/all-geocodes-v2016.xlsx"
tmp<-tempfile()
download.file(url,destfile = tmp,method="curl",mode="wb")
fips<-read_excel(tmp,range="A5:G43939")
fips<-fips[fips$`Area Name (including legal/statistical area description)`=="St. Paul city",]
#Ramsey county is 123, MN is 27
cens<-get_acs(geography = "tract", variables = "B19013_001",state = "MN", county = "Ramsey", geometry = TRUE) #These functions are fantastic!
#geometries are multipolygon sfc_multipolygon
st_crs(cens)
tmp <- tempfile()
download.file("https://information.stpaul.gov/api/geospatial/ykwt-ie3e?method=export&format=Original",tmp,method="curl",mode="wb")
pg2.2<- read_sf(unzip(tmp,"PGrid/PoliceGrid.shp")) #read the original shapefile here, no API token needed.
p<-readOGR(unzip(tmp,"PGrid/PoliceGrid.shp")) #need to identify what is the CRS here.
p<-shapefile(unzip(tmp,"PGrid/PoliceGrid.shp"))
unlink(tmp)
st_crs(pg2.2)
#Police grid map is in WGS84 and I need to convert it to NAD83. CRS = 4326
pg.nad<-st_transform(x = pg2.2, crs = 32610)
#Going to test a map really quick
map.df<-crime2[crime2$incident=="Robbery",] %>% group_by(year,grid) %>% summarize_at("count",sum,na.rm=T)
names(pg2.2)<-c("DIST","grid","geometry")
pg2.3<-left_join(pg2.2,map.df)
pg2.3<-pg2.3[pg2.3$year==2018,c(2:5)]
tmap_mode("view")
tm_shape(pg2.3) + tm_polygons(col="count",n=7,group="grid",
alpha=0.45, title="Robbery Incidents by Grid 2018"
)
##leaflet(cens) %>%
# addPolygons(group="GEOID") %>%
# addTiles(group = "OSM") %>%
# addProviderTiles("CartoDB.Voyager")