library(sf)
library(areal)
library(dplyr)
= c("WHITE", "BLACK", "ASIAN", "HISPANIC", "AM", "OTHER") list_race
Celem ćwiczenia jest przeanalizowanie jak zmieniały się typy struktury rasowo-etnicznej ludności w wybranym hrabstwie pomiędzy dwoma wybranymi latami. W ćwiczeniu zostaną wykorzystane dane zagregowane do obszarów spisowych (ang. census tract) dla District of Columbia dla roku 1990 oraz 2020.
Ćwiczenie składa się z kilku etapów:
Granice obszarów spisowych zmieniały się pomiędzy latami. Aby analiza porównawcza była możliwa, dane trzeba najpierw prztransformować do tego samego zestawu granic. W nieniejszym ćwiczeniu dane z roku 1990 zostaną przetransformowane do granic obwodów spisowych z 2020 roku. W tym celu zostanie wykorzystana metoda powierzchniowo-wagowa.
Używając danych na poziomie obszarów spisowych, proszę przeanalizować zmiany typów zróżnicowania rasowego w wybranym hrabstwie, pomiędzy latami 1990-2000, 2000-2010, 2010-2020, 1990-2020. Wyniki proszę przedstawić w postaci map oraz tabel.
Uwaga! Dane z lat 1990, 2000, 2010 należy przetransformować do granic z 2020 roku. Wykonanie tego ćwiczenia wymaga zestawienia danych wszystkich członków grupy. Dane wykorzystane w tym ćwiczeniu zostały utworzone w trakcie ćwiczenia 2 i zapisane w pliku *[nazwa_hrabstwa]_attr.gpkg* w folderze data/dane_geo.
Uwaga! W przypadku grup 3 osobowych analizujących dane z lat 1990, 2000, 2010 - dane z lat 1990 i 2000 należy przetransformować do granic z 2010 roku. Analizę porównawczą należy wykonać dla lat 1990-2000, 2000-2010, 1990-2010.
Uwaga! W przypadku grup 3 osobowych analizujących dane z lat 1990, 2000, 2010, analizę porównawczą należy wykonać dla lat 1990-2000, 2000-2010, 1990-2010.
Proszę wykonać raport zawierający:
mapy dla analizowanych lat pokazujące rozmieszczenie typów struktury rasowo-etnicznej w hrabstwie.
mapy zmian pokazujące zmiany między latami 1990-2000, 2000-2010, 2010-2020, 1990-2020.
wyniki numeryczne:
komentarz zamieszczonych wyników (ok. 0,5 strony tekstu).
Standaryzacja danych między latami
Poniższy przykład (część 3) pokazuje, w jaki sposób wykonać
transformację danych wykorzystując pakiet areal
w R. W
przykładzie dane z 1990 roku zostały przetransformowane do granic z roku
2020.
Klasyfikację dla lat 1990 i 2020 należy wykonać samodzielnie na podstawie zamieszczonych poniżej kryteriów. Wyniki klasyfikacji należy zapisać jako plik ESRI Shapefile zawierający następujące kolumny: GISJOIN - identyfikator obszaru spisowego, cls90 - klasyfikację dla roku 1990, cls20 - klasyfikację dla roku 2020.
Poniżej (część 5) pokazano przykład, w jaki sposób w R obliczyć liczbę obszarów według typów struktury rasowo-etnicznej oraz w jaki sposób obliczyć macierz przejść.
Uwaga! W pierwszym kroku należy połączyć dane atrybutowe oraz dane przestrzenne na poziomie obszarów spisowych (ang. census_tract). Etap łączenia został wykonany w ćwiczeniu 2. Używając zatem danych zapisanych w pliku *[nazwa_hrabstwa]_attr.gpkg* w folderze data/dane_geo można ten etap pominąć.
= read.csv("data/dane_attr/dc_tract_1990_attr.csv")
csv90 = read.csv("data/dane_attr/dc_tract_2020_attr.csv") csv20
= st_read("data/dane_geo/dc.gpkg", layer = 'tract_1990') fshp90
## Reading layer `tract_1990' from data source
## `/home/anna/DYDAKTYKA/Analiza_geoinformacyjna/cwiczenia2022_23/na_www/data/dane_geo/dc.gpkg'
## using driver `GPKG'
## Simple feature collection with 192 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 1610795 ymin: 1915125 xmax: 1629412 ymax: 1936100
## CRS: 5070
= st_read("data/dane_geo/dc.gpkg", layer = 'tract_2020') fshp20
## Reading layer `tract_2020' from data source
## `/home/anna/DYDAKTYKA/Analiza_geoinformacyjna/cwiczenia2022_23/na_www/data/dane_geo/dc.gpkg'
## using driver `GPKG'
## Simple feature collection with 206 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 1610830 ymin: 1915291 xmax: 1629412 ymax: 1936182
## CRS: 5070
= merge(fshp90[, "GISJOIN"], csv90[,c("GISJOIN_T", list_race, "POP")], by.x = "GISJOIN", by.y = "GISJOIN_T")
dat90 = merge(fshp20[, "GISJOIN"], csv20[,c("GISJOIN_T", list_race, "POP")], by.x = "GISJOIN", by.y = "GISJOIN_T") dat20
st_write(dat90, "data/dane_geo/dc_attr.gpkg", layer = 'tract_1990', append = TRUE)
## Updating layer `tract_1990' to data source `data/dane_geo/dc_attr.gpkg' using driver `GPKG'
## Updating existing layer tract_1990
## Writing 191 features with 8 fields and geometry type Multi Polygon.
st_write(dat20, "data/dane_geo/dc_attr.gpkg", layer = 'tract_2020', append = TRUE)
## Updating layer `tract_2020' to data source `data/dane_geo/dc_attr.gpkg' using driver `GPKG'
## Updating existing layer tract_2020
## Writing 206 features with 8 fields and geometry type Multi Polygon.
Metoda powierzchniowo-wagowa (ang. areal interpolation, areal weightening interpolation) to rodzaj interpolacji powierzchniowej używanej do przekształcania danych geograficznych z jednego zestawu granic w inny. W metodzie tej wartości cechy (np. liczby ludności) przypisywane są do jednostki docelowej proporcjonalnie do odsetka powierzchni, w jakiej jednostka źródłowa pokrywa jednostkę docelową. Jednostką źródłową (ang. source zone) jest pierwotny zestaw granic, w jakim zagregowana jest dana cecha. Jednostką docelową (ang. target zone) jest wynikowy zestaw granic, do którego rozkładamy daną cechę (może nią być inny zestaw granic poligonowych lub siatka komórek rastra).
Metoda powierzchniowo-wagowa została zaimplementowana w pakiecie
areal
w R. Dokumentacja pakietu dostępna jest na stronie:
https://slu-opengis.github.io/areal/.
Obliczenia z wykorzystaniem metody powierzchniowo-wagowej można
wykonać także w oprogramowaniu GIS lub za pomocą funkcji przestrzennych
z pakietu sf. Jednakże, pakiet areal
ma dodatkowo wbudowane
funkcje sprawdzania (walidacji) wyników oraz radzenia sobie w sytuacji,
gdy granice idealnie się nie pokrywają.
Przykład jak wykonać transformację danych między dwoma zestawami granic został dokładnie opisany w dokumentacji pakietu: https://slu-opengis.github.io/areal/articles/areal-weighted-interpolation.html
W niniejszym ćwiczeniu dane dotyczące struktury rasowo-etnicznej z roku 1990 zostały przetransformowane do granic obszarów spisowych z roku 2020.
#target units to granice z 2020 roku
<- dat20[, c("GISJOIN")]
target_units
#source units to dane w granicach z 1990 roku
<- dat90[, c("GISJOIN", list_race, "POP")]
source_units #Uwaga po zapisaniu do GPKG kolumna z geometrią nazywa się geom, a nie geometry. Zatem używając danych wczytanych z pliku gpkg, poniższa linijka powinna mieć zatem postać colnames(source_units) <- c("GJOIN", list_race, "POP", "geom")
colnames(source_units) <- c("GJOIN", list_race, "POP", "geometry")
head(target_units)
## Simple feature collection with 6 features and 1 field
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 1614193 ymin: 1924396 xmax: 1617196 ymax: 1928676
## CRS: 5070
## GISJOIN geometry
## 1 G1100010000101 MULTIPOLYGON (((1617191 192...
## 2 G1100010000102 MULTIPOLYGON (((1616328 192...
## 3 G1100010000201 MULTIPOLYGON (((1614894 192...
## 4 G1100010000202 MULTIPOLYGON (((1615456 192...
## 5 G1100010000300 MULTIPOLYGON (((1615387 192...
## 6 G1100010000400 MULTIPOLYGON (((1615334 192...
head(source_units)
## Simple feature collection with 6 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 1614217 ymin: 1924368 xmax: 1617260 ymax: 1929158
## CRS: 5070
## GJOIN WHITE BLACK ASIAN HISPANIC AM OTHER POP
## 1 G11000100001 4175 147 145 233 3 3 4706
## 2 G1100010000201 2860 371 315 266 7 17 3836
## 3 G1100010000202 3267 60 121 181 6 4 3639
## 4 G1100010000310 5046 186 208 284 8 6 5738
## 5 G11000100004 1073 48 78 114 5 1 1319
## 6 G1100010000501 2332 143 107 183 1 8 2774
## geometry
## 1 MULTIPOLYGON (((1616746 192...
## 2 MULTIPOLYGON (((1614942 192...
## 3 MULTIPOLYGON (((1615496 192...
## 4 MULTIPOLYGON (((1615396 192...
## 5 MULTIPOLYGON (((1615383 192...
## 6 MULTIPOLYGON (((1616620 192...
Metoda powierzchniowo-wagowa składa się z kilku etapów. W pakiecie
areal
można wykonać obliczenia krok po kroku (tak jak zostało to opisane w etapach 1 - 4, lub też użyć funkcjiaw_interpolate()
, która umożliwia wykonanie wszystkich etapów obliczeń od razu dla kilku zmiennych.
Aby wykonać standaryzację danych można pominąć etap 1 - etap 4 oraz przejść od razu do punktu 3.1.1 Funkcja aw_interpolate()
Funkcja aw_intersect()
zostanie wykorzystana do
przecięciea granic jednostek źródłowych (source_units) oraz jednostek
docelowych (target_units). Funkcja ta automatycznie oblicza powierzchnię
każdego poligonu powstałego z przecięcia.
<- aw_intersect(target_units,
intersect source = source_units,
areaVar = "area")
head(intersect)
## Simple feature collection with 6 features and 10 fields
## Geometry type: GEOMETRY
## Dimension: XY
## Bounding box: xmin: 1615456 ymin: 1924399 xmax: 1617196 ymax: 1927355
## CRS: 5070
## GJOIN WHITE BLACK ASIAN HISPANIC AM OTHER POP GISJOIN
## 1 G11000100001 4175 147 145 233 3 3 4706 G1100010000101
## 62 G11000100041 2318 119 142 208 6 6 2799 G1100010000101
## 1.1 G11000100001 4175 147 145 233 3 3 4706 G1100010000102
## 3 G1100010000202 3267 60 121 181 6 4 3639 G1100010000102
## 5 G11000100004 1073 48 78 114 5 1 1319 G1100010000102
## 62.1 G11000100041 2318 119 142 208 6 6 2799 G1100010000102
## area geometry
## 1 203030.190 POLYGON ((1617054 1926648, ...
## 62 2007.179 POLYGON ((1617194 1926546, ...
## 1.1 1606672.169 MULTIPOLYGON (((1616746 192...
## 3 19897.722 MULTIPOLYGON (((1615618 192...
## 5 21968.964 POLYGON ((1616323 1927289, ...
## 62.1 18957.746 POLYGON ((1617054 1926648, ...
W wyniku działania funkcji otrzymamy obiekt zawierający zmienne z dwóch warstw (źródłowej oraz docelowej). Dodana zostanie także kolumna area zawierająca powierzchnię poszczególnych poligonów powstałych z przecięcia.
W drugim etapie obliczane są wagi dla każdego poligonu powstałego z przecięcia granic jednostek źródłowych oraz jednostek docelowych. Wagi dla każdego przecięcia obliczane są według wzoru:
\[W_i = \frac{A_i}{A_j}\]
gdzie
Funkcja aw_total()
obliczy nowe pole zawierające ogólną
powierzchnię wg id jednostek źródłowych. Proszę sprawdzić w dokumentacji
co oznacza argument type=“extensive” oraz
weight=“sum”
<- aw_total(intersect,
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")
head(intersect)
## Simple feature collection with 6 features and 11 fields
## Geometry type: GEOMETRY
## Dimension: XY
## Bounding box: xmin: 1615456 ymin: 1924399 xmax: 1617196 ymax: 1927355
## CRS: 5070
## GJOIN WHITE BLACK ASIAN HISPANIC AM OTHER POP GISJOIN
## 1 G11000100001 4175 147 145 233 3 3 4706 G1100010000101
## 2 G11000100041 2318 119 142 208 6 6 2799 G1100010000101
## 3 G11000100001 4175 147 145 233 3 3 4706 G1100010000102
## 4 G1100010000202 3267 60 121 181 6 4 3639 G1100010000102
## 5 G11000100004 1073 48 78 114 5 1 1319 G1100010000102
## 6 G11000100041 2318 119 142 208 6 6 2799 G1100010000102
## area totalArea geometry
## 1 203030.190 1862902.6 POLYGON ((1617054 1926648, ...
## 2 2007.179 781279.2 POLYGON ((1617194 1926546, ...
## 3 1606672.169 1862902.6 MULTIPOLYGON (((1616746 192...
## 4 19897.722 780479.5 MULTIPOLYGON (((1615618 192...
## 5 21968.964 1620228.2 POLYGON ((1616323 1927289, ...
## 6 18957.746 781279.2 POLYGON ((1617054 1926648, ...
Wagi zostaną obliczone używając funkcji aw_weight()
.
<- aw_weight(intersect,
intersect areaVar = "area", #Ai - obliczane przez aw_intersect
totalVar = "totalArea", #Aj, obliczane przez aw_total
areaWeight = "areaWeight") #Wi
head(intersect)
## Simple feature collection with 6 features and 12 fields
## Geometry type: GEOMETRY
## Dimension: XY
## Bounding box: xmin: 1615456 ymin: 1924399 xmax: 1617196 ymax: 1927355
## CRS: 5070
## GJOIN WHITE BLACK ASIAN HISPANIC AM OTHER POP GISJOIN
## 1 G11000100001 4175 147 145 233 3 3 4706 G1100010000101
## 2 G11000100041 2318 119 142 208 6 6 2799 G1100010000101
## 3 G11000100001 4175 147 145 233 3 3 4706 G1100010000102
## 4 G1100010000202 3267 60 121 181 6 4 3639 G1100010000102
## 5 G11000100004 1073 48 78 114 5 1 1319 G1100010000102
## 6 G11000100041 2318 119 142 208 6 6 2799 G1100010000102
## area totalArea geometry areaWeight
## 1 203030.190 1862902.6 POLYGON ((1617054 1926648, ... 0.108985941
## 2 2007.179 781279.2 POLYGON ((1617194 1926546, ... 0.002569093
## 3 1606672.169 1862902.6 MULTIPOLYGON (((1616746 192... 0.862456362
## 4 19897.722 780479.5 MULTIPOLYGON (((1615618 192... 0.025494229
## 5 21968.964 1620228.2 POLYGON ((1616323 1927289, ... 0.013559179
## 6 18957.746 781279.2 POLYGON ((1617054 1926648, ... 0.024265008
\[E_i = V_j \times W_i\] gdzie
Funkcja aw_calculate()
obliczy szacowaną liczbę osób dla
każdego poligonu \(i\) powstałego z
przecięcia (kolumna POPnew). Wartość ta jest obliczana poprzez
wymnożenie wartości POP oraz areaWeight.
#Jesli argument newVar nie zostanie podany, zostanie nadpisana wartosc w kolumnie value.
<- aw_calculate(intersect,
intersect value = "POP", #kolumna do przeliczenia
newVar = "POPnew", #nazwa nowej kolumny z przeliczonymi danymi.
areaWeight = "areaWeight") #wagi obliczone przez aw_weight()
head(intersect)
## Simple feature collection with 6 features and 13 fields
## Geometry type: GEOMETRY
## Dimension: XY
## Bounding box: xmin: 1615456 ymin: 1924399 xmax: 1617196 ymax: 1927355
## CRS: 5070
## GJOIN WHITE BLACK ASIAN HISPANIC AM OTHER POP GISJOIN
## 1 G11000100001 4175 147 145 233 3 3 4706 G1100010000101
## 2 G11000100041 2318 119 142 208 6 6 2799 G1100010000101
## 3 G11000100001 4175 147 145 233 3 3 4706 G1100010000102
## 4 G1100010000202 3267 60 121 181 6 4 3639 G1100010000102
## 5 G11000100004 1073 48 78 114 5 1 1319 G1100010000102
## 6 G11000100041 2318 119 142 208 6 6 2799 G1100010000102
## area totalArea geometry areaWeight POPnew
## 1 203030.190 1862902.6 POLYGON ((1617054 1926648, ... 0.108985941 512.88784
## 2 2007.179 781279.2 POLYGON ((1617194 1926546, ... 0.002569093 7.19089
## 3 1606672.169 1862902.6 MULTIPOLYGON (((1616746 192... 0.862456362 4058.71964
## 4 19897.722 780479.5 MULTIPOLYGON (((1615618 192... 0.025494229 92.77350
## 5 21968.964 1620228.2 POLYGON ((1616323 1927289, ... 0.013559179 17.88456
## 6 18957.746 781279.2 POLYGON ((1617054 1926648, ... 0.024265008 67.91776
\[G_k = \sum E_{ik}\] gdzie
W ostatnim kroku przeliczone dane (kolumna POPnew obliczona
przez funkcje aw_calculate()
) są agregowane na podstawie id
jednostek docelowych.
<- aw_aggregate(intersect,
result target = target_units, #granice docelowe
tid = GISJOIN, #id docelowe (target units)
interVar = "POPnew") #nazwa przeliczonej zmiennej
aw_interpolate()
Wszystkie 4 etapy można wykonać za pomoca funkcji
aw_interpolate()
. Funkcja aw_interpolate()
umożliwia także wykonanie interpolacji dla kilku zmiennych.
<- aw_interpolate(target_units, #granice jednostek docelowych (target_units)
result90_20 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
Obiekt result90_20 zawiera dane z roku 1990 przeliczone do granic z roku 2020.
Obiekt bnd90_20 zawiera granice z 2020 roku z dołączonymi danymi atrybutowymi z 1990 roku przeliczonymi do granic obszarów spisowych z 2020 roku.
<- merge(dat20[, c("GISJOIN")], result90_20, by = "GISJOIN") bnd90_20
Zapisanie danych atrybutowych z 1990 roku w granicach z 2020 roku (przy wykonywaniu analizy dla innego roku podać nazwę layer = tract_areal_[rok])
st_write(bnd90_20, "data/dane_geo/dc_tract_areal.gpkg", layer = "tract_areal_1990", append = TRUE)
## Updating layer `tract_areal_1990' to data source `data/dane_geo/dc_tract_areal.gpkg' using driver `GPKG'
## Updating existing layer tract_areal_1990
## Writing 206 features with 7 fields and geometry type Multi Polygon.
Dodanie danych z 2020 roku
st_write(dat20, "data/dane_geo/dc_tract_areal.gpkg", layer = "tract_areal_2020", append = TRUE)
## Updating layer `tract_areal_2020' to data source `data/dane_geo/dc_tract_areal.gpkg' using driver `GPKG'
## Updating existing layer tract_areal_2020
## Writing 206 features with 8 fields and geometry type Multi Polygon.
Wykorzystując dane dla roku 1990 (przeliczone do granic z roku 2020) oraz dla roku 2020 (oryginalne dane wczytane z pliku - obiekt dat20) należy sklasyfikować dane do 9 typów struktury rasowo-etnicznej. Poszczególne typy należy oznaczyć odpowiednimi skrótami literowymi - WL, WM itd. (skróty podane są w nawiasie).
Uwaga! Poniżej nie podano kodu w jaki sposób wykonać reklasyfikację. Ten etap trzeba wykonać samodzielnie.
Wynik wykonanej klasyfikacji należy zapisać w folderze results, w pliku [nazwa_hrabstwa]_klasyfikacja9typow.gpkg, layer = cls1990_2020 lub inne lata dla których wykonano zadanie.
Kryterium: grupa dominująca stanowi >80% ogólnej liczby ludności.
Kryterium: grupa dominująca stanowi <50%-80%> ogólnej liczby ludności.
Kryterium: żadna z grup nie przekracza 50%.
Obiekt cls_df zawiera wyniki klasyfikacji struktury rasowo-etnicznej w roku 1990 oraz 2020.
<- st_read("results/dc_klasyfikacja9typow.gpkg", layer = "cls1990_2020") cls_df
## Reading layer `cls1990_2020' from data source
## `/home/anna/DYDAKTYKA/Analiza_geoinformacyjna/cwiczenia2022_23/na_www/results/dc_klasyfikacja9typow.gpkg'
## using driver `GPKG'
## Simple feature collection with 206 features and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 1610830 ymin: 1915291 xmax: 1629412 ymax: 1936182
## CRS: 5070
head(cls_df)
## Simple feature collection with 6 features and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 1614193 ymin: 1924396 xmax: 1617196 ymax: 1928676
## CRS: 5070
## GISJOIN cls90 cls20 geom
## 1 G1100010000101 WL WM MULTIPOLYGON (((1617191 192...
## 2 G1100010000102 WL WM MULTIPOLYGON (((1616328 192...
## 3 G1100010000201 WM WM MULTIPOLYGON (((1614894 192...
## 4 G1100010000202 WL WM MULTIPOLYGON (((1615456 192...
## 5 G1100010000300 WL WM MULTIPOLYGON (((1615387 192...
## 6 G1100010000400 WL WM MULTIPOLYGON (((1615334 192...
table(cls_df$cls90)
##
## AM BL BM HD WL WM
## 1 110 25 15 32 23
Funkcja prob.table()
pozwala na obliczenie odsetka
obszarów według danego typu, jeśli wynik wymnożymy przez 100 otrzymamy
procenty. Funkcja round()
zaokrągli wynik do określonej
liczby miejsc po przecinku
<- prop.table(table(cls_df$cls90))*100
tab90 round(tab90, 1)
Ten sam wynik co powyżej (udział procentowy) otrzymamy także wykonując następujące obliczenia:
table(cls_df$cls90)/sum(table(cls_df$cls90)))*100 (
Wyniki w dwóch powyższych tabelach pokazują, że w roku 1990 75 obszarów (co stanowi 38.5%) zostało sklasyfikowanych jako typ o dużym zróżnicowaniu (żadna grupa nie stanowiła powyżej 50%), a 69 obszarów (35.4%) jako obszary średnio zróżnicowane zdominowane przez białych.
table(cls_df$cls20)
##
## BL BM HD WL WM
## 48 29 43 1 85
<- prop.table(table(cls_df$cls20))*100
tab20 round(tab20, 1)
##
## BL BM HD WL WM
## 23.3 14.1 20.9 0.5 41.3
<- table(cls_df$cls90, cls_df$cls20)
tab tab
##
## BL BM HD WL WM
## AM 0 0 0 0 1
## BL 48 23 26 0 13
## BM 0 4 13 0 8
## HD 0 0 4 0 11
## WL 0 0 0 1 31
## WM 0 2 0 0 21
W powyższej tabeli wiersz zawiera dane dla roku 1990, a kolumny dla roku 2020. Np. wartość 12 (ostatni rząd oraz 4 kolumna) oznacza, że 12 obszarów spisowych sklasyfikowanych w roku 1990 jako WM (mało zróżnicowane zdominowane przez białych) w roku 2020 zmieniło swój typ na HD (duże zróżnicowanie, żadna grupa nie stanowi więcej niż 50%).
Jeśli zsumujemy wartości w wierszach (rowSums) to otrzymamy liczbę
obszarów spisowych według typów w 1990 roku (to samo otrzymamaliśmy
wykonując polecenie table(cls_df$cls90)
)
rowSums(tab)
## AM BL BM HD WL WM
## 1 110 25 15 32 23
Jeśli zsumujemy wartości w wierszach (colSums) to otrzymamy liczbę
obszarów spisowych według typów w 2020 roku (to samo otrzymamaliśmy
wykonując polecenie table(cls_df$cls10)
)
colSums(tab)
## BL BM HD WL WM
## 48 29 43 1 85
Poniższa tabela przedstawia odsetek obszarów spisowych (np. 16.4% w 1990 zostało sklasyfikowanych jako HD, a w 2020 jako AM; 23.6% obszarów spisowych w roku 1990 oraz 20120 było sklasyfikowanych jako WM).
round(prop.table(table(cls_df$cls90, cls_df$cls20))*100, 1)
##
## BL BM HD WL WM
## AM 0.0 0.0 0.0 0.0 0.5
## BL 23.3 11.2 12.6 0.0 6.3
## BM 0.0 1.9 6.3 0.0 3.9
## HD 0.0 0.0 1.9 0.0 5.3
## WL 0.0 0.0 0.0 0.5 15.0
## WM 0.0 1.0 0.0 0.0 10.2
Dwie pierwsze ryciny pokazują przestrzenne rozmieszczenie typów struktury rasowo-etnicznej z zastosowaniem domyślnych kolorów (ten sam typ ma różny kolor w roku 1990 oraz w roku 2020). Mapy te nie nadają się do porównania między latami. Aby mapy były porównywalne muszą być wyświetlane w tych samych kolorach (dany typ musi mieć ten sam kolor na obu mapach)
plot(cls_df["cls90"])
plot(cls_df["cls20"])
Pakiet colorpicker
dostarcza narzędzie do wyboru kolorów
(https://cran.r-project.org/web/packages/colourpicker/vignettes/colourpicker.html)
Po zainstalowaniu pakietu należy wybrać Addins - Colour Picker.
Wektor cls_color zawiera kolory przypisane poszczególnym typom struktury rasowo-etnicznej (Uwaga! Typy muszą być wymienione alfabetycznie).
<- c("AL"= "#CD5555", "AM"= "#FF6A6A", "BL"= "#006400", "BM"= "#32CD32", "HD"= "#8F8F8F", "HL"= "#5D478B", "HM"= "#9370DB", "WL"= "#FF8C00", "WM"= "#FFD700") cls_color
W następnym kroku należy wybrać te typy struktury, które wystąpiły w danym roku.
<- cls_color[names(cls_color)%in%unique(cls_df$cls90)]
col90 <- cls_color[names(cls_color)%in%unique(cls_df$cls20)] col20
Poniższy przykład pokazuje, jak w R wyświetlić dane dla 1990 oraz 2020 roku w tej samej palecie.
library(ggplot2)
library(patchwork)
<- ggplot(data = cls_df) +
p1 geom_sf(aes(fill = cls90)) +
scale_fill_manual(values = col90) +
labs(title = "District of Columbia, 1990") +
theme_bw() +
theme(legend.position="bottom")
<- ggplot(data = cls_df) +
p2 geom_sf(aes(fill = cls20)) +
scale_fill_manual(values = col20) +
labs(title = "District of Columbia, 2020") +
theme_bw() +
theme(legend.position="bottom")
#wyswietlenie wykresow obok siebie
+ p2 p1
Ten etap zostanie wykonany w oprogramowaniu QGIS.
Powyższy przykład prezentował dwie mapy pokazujące typy struktury rasowo-etnicznej - jedną dla roku 1990 oraz jedną dla 2020 roku.
Poniżej zmiany typów struktury rasowo-etnicznej zostaną przedstawione na jednej mapie:
obszary wypełnione jednolitym kolorem nie zmieniły się.
obszary wypełnione szrafem zmieniły się pomiędzy latami:
W QGIS należy wykonać stylizację (patrz rycina poniżej).
Etap 1: Wczytać do QGIS warstwę z klasyfikacją. Przypisać kolory do 9 klas na podstawie zmiennej cls90 zawierającej klasyfikację dla 1990 roku (można wykorzystać plik zmiany.qml).
Rycina przedstawia typy zróżnicowania rasowo-etnicznego w 1990 roku.
Etap 2: - Skopiować warstwę. Dla górnej warstwy zmienić klasę na cls20.
Etap 3: Zmiana stylizacji warstwy dla 2020 roku.
Wybrać Symbol. Zmienić sposób wyświetlania z Simple Fill do Line Pattern Fill. Zmienić odległość między liniami (Spacing) na 3. Kliknąć na Simple Line i wpisać grubość linii 1.
W Wyniku otrzymamy mapę: