cw6

mail@pastecode.io avatar
unknown
r
3 years ago
16 kB
5
Indexable
Never
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