Untitled
unknown
r
4 years ago
8.2 kB
9
Indexable
library(sf) library(areal) library(dplyr) setwd("cw4") list_race = c("WHITE", "BLACK", "ASIAN", "HISPANIC", "INDIAN", "OTHER_RACE") dat90 = st_read("data/san_francisco_tracts_1990_attr.shp") dat00 = st_read("data/san_francisco_tracts_2000_attr.shp") dat10 = st_read("data/san_francisco_tracts_2010_attr.shp") target_units <- dat10[, c("GISJOIN")] source_units <- dat90[, c("GISJOIN", list_race, "POP")] source_units2 <- dat00[, c("GISJOIN", list_race, "POP")] colnames(source_units) <- c("GJOIN", list_race, "POP", "geometry") colnames(source_units2) <- c("GJOIN", list_race, "POP", "geometry") intersect <- aw_intersect(target_units, source = source_units, areaVar = "area") intersect <- aw_total(intersect, source = source_units, #granice zrodlowe id = GJOIN, #id jednostek zrodlowych areaVar = "area", #nazwa pola z powierzchnia jednostek z przeciecia, Ai totalVar = "totalArea", type = "extensive", weight = "sum") intersect <- aw_weight(intersect, areaVar = "area", #Ai - obliczane przez aw_intersect totalVar = "totalArea", #Aj, obliczane przez aw_total areaWeight = "areaWeight") #Wi intersect <- aw_calculate(intersect, value = "POP", #kolumna do przeliczenia newVar = "POPnew", #nazwa nowej kolumny z przeliczonymi danymi. areaWeight = "areaWeight") result <- aw_aggregate(intersect, target = target_units, #granice docelowe tid = GISJOIN, #id docelowe (target units) interVar = "POPnew") #nazwa przeliczonej zmiennej result90_10 <- aw_interpolate(target_units, #granice jednostek docelowych (target_units) tid = GISJOIN, #id jednostek docelowych source = source_units, #dane zrodlowe sid = GJOIN, #id jednostek zrodlowych (source units) weight = "sum", output = "tibble", #wynik w jako tabela, moze tez być obiekt sf extensive = list_race) #lista zmiennych do przeliczenia result00_10 <- aw_interpolate(target_units, #granice jednostek docelowych (target_units) tid = GISJOIN, #id jednostek docelowych source = source_units2, #dane zrodlowe sid = GJOIN, #id jednostek zrodlowych (source units) weight = "sum", output = "tibble", #wynik w jako tabela, moze tez być obiekt sf extensive = list_race) #lista zmiennych do przeliczenia write.csv(result90_10, "result90_10.csv", row.names = FALSE) write.csv(result00_10, "result00_10.csv", row.names = FALSE) dat10 cls_df <- st_read("data/klasyfikacja_1990_2010.shp") result90_10$ASIAN = as.integer(result90_10$ASIAN) result90_10$BLACK = as.integer(result90_10$BLACK) result90_10$HISPANIC = as.integer(result90_10$HISPANIC) result90_10$INDIAN = as.integer(result90_10$INDIAN) result90_10$OTHER_RACE = as.integer(result90_10$OTHER_RACE) result90_10$WHITE = as.integer(result90_10$WHITE) result90_10$POP = result90_10$ASIAN + result90_10$BLACK + result90_10$HISPANIC + result90_10$INDIAN+result90_10$OTHER_RACE+result90_10$WHITE View(result90_10) result00_10$ASIAN = as.integer(result00_10$ASIAN) result00_10$BLACK = as.integer(result00_10$BLACK) result00_10$HISPANIC = as.integer(result00_10$HISPANIC) result00_10$INDIAN = as.integer(result00_10$INDIAN) result00_10$OTHER_RACE = as.integer(result00_10$OTHER_RACE) result00_10$WHITE = as.integer(result00_10$WHITE) result00_10$POP = result00_10$ASIAN + result00_10$BLACK + result00_10$HISPANIC + result00_10$INDIAN+result00_10$OTHER_RACE+result00_10$WHITE View(result00_10) res9_1 = result90_10 %>% summarize(GISJOIN, ASIAN = ASIAN/POP, BLACK = BLACK/POP, HISPANIC = HISPANIC/POP, INDIAN = INDIAN/POP, OTHER_RACE = OTHER_RACE/POP, WHITE = WHITE/POP) res0_1 = result00_10 %>% summarize(GISJOIN, ASIAN = ASIAN/POP, BLACK = BLACK/POP, HISPANIC = HISPANIC/POP, INDIAN = INDIAN/POP, OTHER_RACE = OTHER_RACE/POP, WHITE = WHITE/POP) View(res9_1) View(res0_1) res_9_1 = res9_1 %>% mutate(cls90 = case_when( ASIAN > 0.8 ~ "AL", BLACK > 0.8 ~ "BL", HISPANIC > 0.8 ~ "HL", WHITE > 0.8 ~ "WL", WHITE >= 0.5 & WHITE <= 0.8 ~ "WM", BLACK >= 0.5 & BLACK <= 0.8 ~ "BM", ASIAN >= 0.5 & ASIAN <= 0.8 ~ "AM", HISPANIC >= 0.5 & HISPANIC <= 0.8 ~ "HM", HISPANIC < 0.5 | ASIAN < 0.5 | BLACK < 0.5 | WHITE < 0.5 ~ "HD" )) res_0_1 = res0_1 %>% mutate(cls00 = case_when( ASIAN > 0.8 ~ "AL", BLACK > 0.8 ~ "BL", HISPANIC > 0.8 ~ "HL", WHITE > 0.8 ~ "WL", WHITE >= 0.5 & WHITE <= 0.8 ~ "WM", BLACK >= 0.5 & BLACK <= 0.8 ~ "BM", ASIAN >= 0.5 & ASIAN <= 0.8 ~ "AM", HISPANIC >= 0.5 & HISPANIC <= 0.8 ~ "HM", HISPANIC < 0.5 | ASIAN < 0.5 | BLACK < 0.5 | WHITE < 0.5 ~ "HD" )) cls_df <- st_read("data/klasyfikacja_1990_2010.shp") cls_df$cls90 = res_9_1$cls90 cls_df$cls00 = res_0_1$cls00 View(cls_df) table(cls_df$cls90) tab90 <- prop.table(table(cls_df$cls90))*100 a = round(tab90, 1) a = data.frame(rbind(a)) write.csv(a, "udzial_procentowy_90.csv", row.names = FALSE) table(cls_df$cls00) tab90 <- prop.table(table(cls_df$cls00))*100 a = round(tab90, 1) a = data.frame(rbind(a)) write.csv(a, "udzial_procentowy_00.csv", row.names = FALSE) table(cls_df$cls10) tab90 <- prop.table(table(cls_df$cls10))*100 a = round(tab90, 1) a = data.frame(rbind(a)) write.csv(a, "udzial_procentowy_10.csv", row.names = FALSE) dat10 result10_10$ASIAN = as.integer(result00_10$ASIAN) result10_10$BLACK = as.integer(result00_10$BLACK) result10_10$HISPANIC = as.integer(result00_10$HISPANIC) result10_10$INDIAN = as.integer(result00_10$INDIAN) result10_10$OTHER_RACE = as.integer(result00_10$OTHER_RACE) result10_10$WHITE = as.integer(result00_10$WHITE) result10_10$POP = result00_10$ASIAN + result00_10$BLACK + result00_10$HISPANIC + result00_10$INDIAN+result00_10$OTHER_RACE+result00_10$WHITE View(result00_10) res1_1 = dat10 %>% summarize(GISJOIN, ASIAN = ASIAN/POP, BLACK = BLACK/POP, HISPANIC = HISPANIC/POP, INDIAN = INDIAN/POP, OTHER_RACE = OTHER_RACE/POP, WHITE = WHITE/POP) res1_1 = res1_1 %>% mutate(cls10 = case_when( ASIAN > 0.8 ~ "AL", BLACK > 0.8 ~ "BL", HISPANIC > 0.8 ~ "HL", WHITE > 0.8 ~ "WL", WHITE >= 0.5 & WHITE <= 0.8 ~ "WM", BLACK >= 0.5 & BLACK <= 0.8 ~ "BM", ASIAN >= 0.5 & ASIAN <= 0.8 ~ "AM", HISPANIC >= 0.5 & HISPANIC <= 0.8 ~ "HM", HISPANIC < 0.5 | ASIAN < 0.5 | BLACK < 0.5 | WHITE < 0.5 ~ "HD" )) cls_df$cls10 = res1_1$cls10 View(cls_df) cls_color <- c("AL"= "#CD5555", "AM"= "#FF6A6A", "BL"= "#006400", "BM"= "#32CD32", "HD"= "#8F8F8F", "HL"= "#5D478B", "HM"= "#9370DB", "WL"= "#FF8C00", "WM"= "#FFD700") col90 <- cls_color[names(cls_color)%in%unique(cls_df$cls90)] col00 <- cls_color[names(cls_color)%in%unique(cls_df$cls00)] col10 <- cls_color[names(cls_color)%in%unique(cls_df$cls10)] library(ggplot2) library(gridExtra) p1 <- ggplot(data = cls_df) + geom_sf(aes(fill = cls90)) + scale_fill_manual(values = col90) + labs(title = "San Francisco, 1990") + theme_bw() + theme(legend.position="bottom") p2 <- ggplot(data = cls_df) + geom_sf(aes(fill = cls00)) + scale_fill_manual(values = col10) + labs(title = "San Francisco, 2000") + theme_bw() + theme(legend.position="bottom") p3 <- ggplot(data = cls_df) + geom_sf(aes(fill = cls10)) + scale_fill_manual(values = col10) + labs(title = "San Francisco, 2010") + theme_bw() + theme(legend.position="bottom") grid.arrange(p1, p2, p3, nrow = 2) png("strukturyrasowe.png") grid.arrange(p1, p2, p3, nrow = 2) dev.off() cls_df write_sf(cls_df, "classification.shp")
Editor is loading...