R - skazani na działce() - Farbowanie Shapefile wielokątów w oparciu o wartości szczeliny

Question

Mam shapefile przedstawiający oddalonych obszarów w Australii, uzyskany z Australian Bureau of Statistics:

http://www.abs.gov.au/AUSSTATS/abs@.nsf/DetailsPage/1270.0.55.005July%202011?OpenDocument

W tym samym adresem URL jest PDF "ASGS Remoteness Structure Edition 2011 Mapy PDF" - Próbuję odtworzyć pierwszą mapę z tego dokumentu PDF.

Czytałem w shapefile i dodaje informacje o kolorze do gniazda data:

ra <- readShapeSpatial("RA_2011_AUST", delete_null_obj = TRUE) 
ra@data$COLOUR <- "#FFFFFF" 
ra@data$COLOUR[(as.numeric(as.character(ra@data$RA_CODE11)) %% 10) == 0] <- "#006837" 
ra@data$COLOUR[(as.numeric(as.character(ra@data$RA_CODE11)) %% 10) == 1] <- "#31A354" 
ra@data$COLOUR[(as.numeric(as.character(ra@data$RA_CODE11)) %% 10) == 2] <- "#78C679" 
ra@data$COLOUR[(as.numeric(as.character(ra@data$RA_CODE11)) %% 10) == 3] <- "#C2E699" 
ra@data$COLOUR[(as.numeric(as.character(ra@data$RA_CODE11)) %% 10) == 4] <- "#FFFFCC" 

Jedyne co pozostaje mi zrobić, to działka mapa! To tutaj utknęłam ...

ra@polygons to lista 35 wielokątów, z których każdy ma szczelinę ID, która jest indeksem ramki danych ra@data. Więc wszystko, co muszę zrobić, to powiedzieć plot(), aby znaleźć kolor w ra@data$COLOUR[ID]. Cóż, niezupełnie. Każdy z 35 wielokątów (klasa "Wielokąty") ma własną listę wielokątów (klasa "Wielokąt"); w sumie jest 6902 wielokątów !!!

Moje zrozumienie plot() polega na tym, że muszę przekazać mu wektor kolorów w tej samej kolejności, w jakiej będą drukowane wielokąty. Dlatego uważam, że będę musiał stworzyć wektor o długości 6902 z każdym elementem utrzymującym wartość koloru dla powiązanego wielokąta. Jak ja to robię?

Byłoby to łatwe, gdyby wielokąty zostały naniesione w porządku, ale tak nie jest. Każdy z 35 wielokątów ma szczelinę plotOrder, która jest wektorem całkowitym, więc wektor koloru będzie prawdopodobnie musiał być uporządkowany według wartości w każdym z tych wektorów.

W tym momencie wszystko wydaje się nieco zbyt skomplikowane. Czy jestem całkowicie nie na miejscu?

Dzięki za porady!

Answer 1

Zrobiłeś już całą pracę!

plot(ra, col=ra@data$COLOUR) 

Albo nawet, jak @Spacedman zasugerował:

plot(ra, col=ra$COLOUR) 

I to jest to!

A jeśli chcesz się pozbyć granicach wielokąta:

plot(ra, col=ra$COLOUR, border=NA) 

Edit: próba wyjaśnienia tego zachowania:

Według ?SpatialPointsDataFrame:

SpatialPolygonsDataFrame z domyślnym dopasowaniem identyfikatorów sprawdza nazwy wierszy ramek danych w polach Identyfikatory wielokątów. Muszą się ze sobą zgodzić i być unikatowe (żadne obiekty Wieloboków nie mogą udostępniać identyfikatorów); wiersze ramki danych zostaną ponownie zamówione, jeśli będą potrzebne, aby dopasować identyfikatory wielokątów.

Oznacza to, że wielokąty są uporządkowane zgodnie z ich identyfikatorem, a zatem są w kolejności rzędów ramki danych w gnieździe @data.

Teraz, jeśli spojrzeć na funkcję plot.SpatialPolygons (używając getAnywhere(plot.SpatialPolygons)) są te linie w pewnym momencie:

... 
    polys <- slot(x, "polygons") 
    pO <- slot(x, "plotOrder") 
    if (!is.null(density)) { 
     if (missing(col)) 
      col <- par("fg") 
     if (length(col) != n) 
      col <- rep(col, n, n) 
     if (length(density) != n) 
      density <- rep(density, n, n) 
     if (length(angle) != n) 
      angle <- rep(angle, n, n) 
     for (j in pO) .polygonRingHoles(polys[[j]], border = border[j], 
      xpd = xpd, density = density[j], angle = angle[j], 
      col = col[j], pbg = pbg, lty = lty, ...) 
    } 
... 

Wektor wprowadzane do col ma takiej samej kolejności jak gniazda polygons i stąd jako ID. plotOrder służy do indeksowania wszystkich z nich w ten sam sposób.