cw6
unknown
r
4 years ago
16 kB
9
Indexable
install.packages("raceland") library(raster) library(dplyr) library(raceland) library(sf) library(raster) race_colors = c("#F16667", "#6EBE44", "#7E69AF", "#C77213", "#F8DF1D") biv_colors = c("11" = "#e8e8e8", "12" = "#e4acac", "13" = "#c85a5a", "21" = "#b0d5df", "22" = "#ad9ea5", "23" = "#985356", "31" = "#64acbe", "32"= "#627f8c", "33" = "#574249") ent_breaks = c(seq(0, 2, by = 0.25), log2(nlayers(race_raster))) mut_breaks = seq(0, 1, by = 0.1) setwd("../") getwd() lr <- list.files("cwiczenie_6/dane/cwiczenie", full.names = TRUE) lr library(raster) rast <- stack(lr) rast #zsumowanie dwóch warstw rast[["other"]] = rast[["other"]] + rast[["am"]] #usunięcie niepotrzebnych warstw rast <- dropLayer(rast, "am") rast[["asian"]] = rast[["asian"]] + rast[["pi"]] rast <- dropLayer(rast, "pi") #zdefiniowanie zasiegu do ktorego docinamy warstwe b = as(extent(-2284890, -2270940, 1947390, 1961430), 'SpatialPolygons') #przypisanie ukladu wspolrzednych takiego jak w oryginalnej warstwie crs(b) <- crs(rast) # x - dane rastrowe, n - liczba realizacji real_raster = create_realizations(x = race_raster, n = 100) # wyswietlenie realizacji 1-5. plot(real_raster[[1:5]], col = c("#F16667", "#6EBE44", "#7E69AF", "#C77213", "#F8DF1D")) plot_realization(x = real_raster[[1]], y = race_raster, hex = race_colors) #real raster - krajobraz rasowy, race_raster - obiekt zawierający dane rastrowe, każda komórka ma przypisaną gęstość zaludnienia dla danej rasy. #window_size =10 oznacza, że lokalna gęstość zaludnienia liczona jest dla obszarów 10x10 komórek. dens_raster = create_densities(real_raster, race_raster, window_size = 10) metr_df = calculate_metrics(x = real_raster, #obiekt z realizacjami krajobrazu rasowego w = dens_raster, #obiekt z lokalną gęstością zaludnienia neighbourhood = 4, fun = "mean", size = NULL, threshold = 1) grid_sf = create_grid(real_raster, size = 240) #size 240 oznacza 240x240 komórek. plot_realization(real_raster[[1]], race_raster, hex = race_colors) plot(st_geometry(grid_sf), add = TRUE, lwd = 1) metr_df_1 = calculate_metrics(x = real_raster, w = dens_raster, neighbourhood = 4, fun = "mean", size = 240, # TU NALEZY ZDEFINIOWAC WIELKOSC OCZKA SIATKI - dla jakiej skali wykonujemy obliczenia. threshold = 0.6) metr_df_1[metr_df_1$realization == 1, ] by_row_col <- group_by(metr_df_1, row, col) smr <- summarize(by_row_col, ent_mean = mean(ent, na.rm = TRUE), ent_sd = sd(ent, na.rm = TRUE), mutinf_mean = mean(mutinf, na.rm = TRUE), mutinf_sd = sd(mutinf, na.rm = TRUE) ) smr attr_grid = dplyr::left_join(grid_sf, smr, by = c("row", "col")) plot(attr_grid["ent_mean"], breaks = ent_breaks, key.pos = 1, pal = rev(RColorBrewer::brewer.pal(length(ent_breaks) - 1, name = "RdBu")), # pal = grDevices::hcl.colors(length(ent_breaks) - 1, palette = "Blue-Red"), bty = "n", main = "Racial diversity (Entropy)") plot(attr_grid["mutinf_mean"], breaks = mut_breaks, key.pos = 1, pal = rev(RColorBrewer::brewer.pal(length(mut_breaks) - 1, name = "RdBu")), # pal = grDevices::hcl.colors(length(mut_breaks) - 1, palette = "Blue-Red"), bty = "n", main = "Racial segregation (Mutual information)") bivariate_classification = function(entropy, mutual_information, n) { # calculate bivariate classification nent = log2(n) ent_cat = cut(entropy, breaks = c(0, 0.66, 1.33, nent), labels = c(1, 2, 3), include.lowest = TRUE, right = TRUE) ent_cat = as.integer(as.character(ent_cat)) mut_cat = cut(mutual_information, breaks = c(0, 0.33, 0.66, 1), labels = c(10, 20, 30), include.lowest = TRUE, right = TRUE) mut_cat = as.integer(as.character(mut_cat)) bivar_cls = mut_cat + ent_cat bivar_cls = as.factor(bivar_cls) return(bivar_cls) } attr_grid$bivar_cls = bivariate_classification(entropy = smr$ent_mean, mutual_information = smr$mutinf_mean, n = nlayers(race_raster)) biv_colors = c("11" = "#e8e8e8", "12" = "#e4acac", "13" = "#c85a5a", "21" = "#b0d5df", "22" = "#ad9ea5", "23" = "#985356", "31" = "#64acbe", "32"= "#627f8c", "33" = "#574249") bcat = biv_colors[names(biv_colors)%in%unique(attr_grid$bivar_cls)] plot(attr_grid["bivar_cls"], pal = bcat, main = "Racial diversity and residential segregation") #przyciecie race_raster <- crop(rast, b) ############3 ZADANIA #1 setwd("cwiczenie_6/dane") getwd() file.rename("ca_sanfrancisco_nhw2010myc.tif", "white.tif") file.rename("ca_sanfrancisco_nhb2010myc.tif", "black.tif") file.rename("ca_sanfrancisco_nhas2010myc.tif", "asian.tif") file.rename("ca_sanfrancisco_hispanic2010myc.tif", "hispanic.tif") file.rename("ca_sanfrancisco_nhother2010myc.tif", "other_race.tif") file.rename("ca_sanfrancisco_nhpi2010myc.tif", "pi.tif") file.rename("ca_sanfrancisco_nham2010myc.tif", "am.tif") setwd("../") lr <- list.files("dane/2010",full.names = TRUE) lr library(raster) rast <- stack(lr) rast rast[["other_race"]] = rast[["other_race"]] + rast[["am"]] rast <- dropLayer(rast, "am") rast[["asian"]] = rast[["asian"]] + rast[["pi"]] rast <- dropLayer(rast, "pi") #zdefiniowanie zasiegu do ktorego docinamy warstwe b = as(extent(-2284890, -2270940, 1947390, 1961430), 'SpatialPolygons') #przypisanie ukladu wspolrzednych takiego jak w oryginalnej warstwie crs(b) <- crs(rast) #przyciecie race_raster <- crop(rast, b) #plot(race_raster) # x - dane rastrowe, n - liczba realizacji real_raster = create_realizations(x = race_raster, n = 100) plot(real_raster[[1:5]], col = c("#F16667", "#6EBE44", "#7E69AF", "#C77213", "#F8DF1D")) plot_realization(x = real_raster[[1]], y = race_raster, hex = race_colors) dens_raster = create_densities(real_raster, race_raster, window_size = 10) metr_df = calculate_metrics(x = real_raster, #obiekt z realizacjami krajobrazu rasowego w = dens_raster, #obiekt z lokalną gęstością zaludnienia neighbourhood = 4, fun = "mean", size = NULL, threshold = 1) summarise(metr_df, mean_ent = mean(ent, na.rm = TRUE), sd_ent = sd(ent, na.rm = TRUE), mean_mutinf = mean(mutinf), sd_mutinf = sd(mutinf)) grid_sf = create_grid(real_raster, size = 60) #size 240 oznacza 240x240 komórek. metr_df_1 = calculate_metrics(x = real_raster, w = dens_raster, neighbourhood = 4, fun = "mean", size = 60, # TU NALEZY ZDEFINIOWAC WIELKOSC OCZKA SIATKI - dla jakiej skali wykonujemy obliczenia. threshold = 0.6) plot_realization(real_raster[[1]], race_raster, hex = race_colors) plot(st_geometry(grid_sf), add = TRUE, lwd = 1) metr_df_1[metr_df_1$realization == 1, ] by_row_col <- group_by(metr_df_1, row, col) smr <- summarize(by_row_col, ent_mean = mean(ent, na.rm = TRUE), ent_sd = sd(ent, na.rm = TRUE), mutinf_mean = mean(mutinf, na.rm = TRUE), mutinf_sd = sd(mutinf, na.rm = TRUE) ) smr # to do wyjęcia write.csv(smr, "summarize_60_2010.csv", row.names = FALSE) attr_grid = dplyr::left_join(grid_sf, smr, by = c("row", "col")) # Zróżnicowanie rasowe png("diversity2010_60.jpg") plot(attr_grid["ent_mean"], breaks = ent_breaks, key.pos = 1, pal = rev(RColorBrewer::brewer.pal(length(ent_breaks) - 1, name = "RdBu")), # pal = grDevices::hcl.colors(length(ent_breaks) - 1, palette = "Blue-Red"), bty = "n", main = "Racial diversity (Entropy)") dev.off() # Segregacja rasowa png("segregation2010_60.jpg") plot(attr_grid["mutinf_mean"], breaks = mut_breaks, key.pos = 1, pal = rev(RColorBrewer::brewer.pal(length(mut_breaks) - 1, name = "RdBu")), # pal = grDevices::hcl.colors(length(mut_breaks) - 1, palette = "Blue-Red"), bty = "n", main = "Racial segregation (Mutual information)") dev.off() bivariate_classification = function(entropy, mutual_information, n) { # calculate bivariate classification nent = log2(n) ent_cat = cut(entropy, breaks = c(0, 0.66, 1.33, nent), labels = c(1, 2, 3), include.lowest = TRUE, right = TRUE) ent_cat = as.integer(as.character(ent_cat)) mut_cat = cut(mutual_information, breaks = c(0, 0.33, 0.66, 1), labels = c(10, 20, 30), include.lowest = TRUE, right = TRUE) mut_cat = as.integer(as.character(mut_cat)) bivar_cls = mut_cat + ent_cat bivar_cls = as.factor(bivar_cls) return(bivar_cls) } attr_grid$bivar_cls = bivariate_classification(entropy = smr$ent_mean, mutual_information = smr$mutinf_mean, n = nlayers(race_raster)) biv_colors = c("11" = "#e8e8e8", "12" = "#e4acac", "13" = "#c85a5a", "21" = "#b0d5df", "22" = "#ad9ea5", "23" = "#985356", "31" = "#64acbe", "32"= "#627f8c", "33" = "#574249") bcat = biv_colors[names(biv_colors)%in%unique(attr_grid$bivar_cls)] png("racial2010_60.jpg") plot(attr_grid["bivar_cls"], pal = bcat, main = "Racial diversity and residential segregation") dev.off() # to były typy struktury rasowo etnicznej ################ DLA 120 grid_sf = create_grid(real_raster, size = 120) #size 240 oznacza 240x240 komórek. metr_df_1 = calculate_metrics(x = real_raster, w = dens_raster, neighbourhood = 4, fun = "mean", size = 120, # TU NALEZY ZDEFINIOWAC WIELKOSC OCZKA SIATKI - dla jakiej skali wykonujemy obliczenia. threshold = 0.6) metr_df_1[metr_df_1$realization == 1, ] by_row_col <- group_by(metr_df_1, row, col) smr <- summarize(by_row_col, ent_mean = mean(ent, na.rm = TRUE), ent_sd = sd(ent, na.rm = TRUE), mutinf_mean = mean(mutinf, na.rm = TRUE), mutinf_sd = sd(mutinf, na.rm = TRUE) ) smr # to do wyjęcia write.csv(smr, "summarize_120_2010.csv", row.names = FALSE) attr_grid = dplyr::left_join(grid_sf, smr, by = c("row", "col")) # Zróżnicowanie rasowe png("diversity2010_120.jpg") plot(attr_grid["ent_mean"], breaks = ent_breaks, key.pos = 1, pal = rev(RColorBrewer::brewer.pal(length(ent_breaks) - 1, name = "RdBu")), # pal = grDevices::hcl.colors(length(ent_breaks) - 1, palette = "Blue-Red"), bty = "n", main = "Racial diversity (Entropy)") dev.off() # Segregacja rasowa png("segregation2010_120.jpg") plot(attr_grid["mutinf_mean"], breaks = mut_breaks, key.pos = 1, pal = rev(RColorBrewer::brewer.pal(length(mut_breaks) - 1, name = "RdBu")), # pal = grDevices::hcl.colors(length(mut_breaks) - 1, palette = "Blue-Red"), bty = "n", main = "Racial segregation (Mutual information)") dev.off() bivariate_classification = function(entropy, mutual_information, n) { # calculate bivariate classification nent = log2(n) ent_cat = cut(entropy, breaks = c(0, 0.66, 1.33, nent), labels = c(1, 2, 3), include.lowest = TRUE, right = TRUE) ent_cat = as.integer(as.character(ent_cat)) mut_cat = cut(mutual_information, breaks = c(0, 0.33, 0.66, 1), labels = c(10, 20, 30), include.lowest = TRUE, right = TRUE) mut_cat = as.integer(as.character(mut_cat)) bivar_cls = mut_cat + ent_cat bivar_cls = as.factor(bivar_cls) return(bivar_cls) } attr_grid$bivar_cls = bivariate_classification(entropy = smr$ent_mean, mutual_information = smr$mutinf_mean, n = nlayers(race_raster)) biv_colors = c("11" = "#e8e8e8", "12" = "#e4acac", "13" = "#c85a5a", "21" = "#b0d5df", "22" = "#ad9ea5", "23" = "#985356", "31" = "#64acbe", "32"= "#627f8c", "33" = "#574249") bcat = biv_colors[names(biv_colors)%in%unique(attr_grid$bivar_cls)] png("racial2010_120.jpg") plot(attr_grid["bivar_cls"], pal = bcat, main = "Racial diversity and residential segregation") dev.off() # to były typy struktury rasowo etnicznej ##### ZADANIE 3 ## 1990 getwd() setwd("2000") # 30, 60, 120, 240 getwd() file.rename("ca_sanfrancisco_nhw2000myc.tif", "white.tif") file.rename("ca_sanfrancisco_nhb2000myc.tif", "black.tif") file.rename("ca_sanfrancisco_nhas2000myc.tif", "asian.tif") file.rename("ca_sanfrancisco_hispanic2000myc.tif", "hispanic.tif") file.rename("ca_sanfrancisco_nhother2000myc.tif", "other_race.tif") file.rename("ca_sanfrancisco_nhpi2000myc.tif", "pi.tif") file.rename("ca_sanfrancisco_nham2000myc.tif", "am.tif") setwd("../") lr <- list.files("2000", full.names = TRUE) lr library(raster) rast <- stack(lr) rast rast[["other_race"]] = rast[["other_race"]] + rast[["am"]] rast <- dropLayer(rast, "am") rast[["asian"]] = rast[["asian"]] + rast[["pi"]] rast <- dropLayer(rast, "pi") b = as(extent(-2284890, -2270940, 1947390, 1961430), 'SpatialPolygons') crs(b) <- crs(rast) race_raster <- crop(rast, b) real_raster = create_realizations(x = race_raster, n = 100) plot(real_raster[[1:5]], col = c("#F16667", "#6EBE44", "#7E69AF", "#C77213", "#F8DF1D")) dens_raster = create_densities(real_raster, race_raster, window_size = 10) metr_df = calculate_metrics(x = real_raster, #obiekt z realizacjami krajobrazu rasowego w = dens_raster, #obiekt z lokalną gęstością zaludnienia neighbourhood = 4, fun = "mean", size = NULL, threshold = 1) summarise(metr_df, mean_ent = mean(ent, na.rm = TRUE), sd_ent = sd(ent, na.rm = TRUE), mean_mutinf = mean(mutinf), sd_mutinf = sd(mutinf)) grid_sf_30 = create_grid(real_raster, size = 240) #size 240 oznacza 240x240 komórek. metr_df_1_30 = calculate_metrics(x = real_raster, w = dens_raster, neighbourhood = 4, fun = "mean", size = 240, # TU NALEZY ZDEFINIOWAC WIELKOSC OCZKA SIATKI - dla jakiej skali wykonujemy obliczenia. threshold = 0.6) by_row_col_30 <- group_by(metr_df_1_30, row, col) smr_30 <- summarize(by_row_col_30, ent_mean = mean(ent, na.rm = TRUE), ent_sd = sd(ent, na.rm = TRUE), mutinf_mean = mean(mutinf, na.rm = TRUE), mutinf_sd = sd(mutinf, na.rm = TRUE) ) write.csv(smr_30, "summarize_240_2000.csv", row.names = FALSE) attr_grid = dplyr::left_join(grid_sf_30, smr_30, by = c("row", "col")) png("diversity_2000_240.png") plot(attr_grid["ent_mean"], breaks = ent_breaks, key.pos = 1, pal = rev(RColorBrewer::brewer.pal(length(ent_breaks) - 1, name = "RdBu")), # pal = grDevices::hcl.colors(length(ent_breaks) - 1, palette = "Blue-Red"), bty = "n", main = "Racial diversity (Entropy)") dev.off() png("segregation_2000_240.png") plot(attr_grid["mutinf_mean"], breaks = mut_breaks, key.pos = 1, pal = rev(RColorBrewer::brewer.pal(length(mut_breaks) - 1, name = "RdBu")), # pal = grDevices::hcl.colors(length(mut_breaks) - 1, palette = "Blue-Red"), bty = "n", main = "Racial segregation (Mutual information)") dev.off() attr_grid$bivar_cls = bivariate_classification(entropy = smr_30$ent_mean, mutual_information = smr_30$mutinf_mean, n = nlayers(race_raster)) biv_colors = c("11" = "#e8e8e8", "12" = "#e4acac", "13" = "#c85a5a", "21" = "#b0d5df", "22" = "#ad9ea5", "23" = "#985356", "31" = "#64acbe", "32"= "#627f8c", "33" = "#574249") bcat = biv_colors[names(biv_colors)%in%unique(attr_grid$bivar_cls)] png("racial_2000_240.png") plot(attr_grid["bivar_cls"], pal = bcat, main = "Racial diversity and residential segregation") dev.off() ### 60
Editor is loading...