Usa Rcartogram su un oggetto SpatialPolygonsDataFrame

Question

Sto provando a fare la stessa cosa posta in questa domanda, Cartogram + choropleth map in R, ma partendo da SpatialPolygonsDataFrame e sperando di finire con lo stesso tipo di oggetto.

È possibile salvare l'oggetto come file di forma, utilizzare scapetoad, riaprirlo e riconvertire, ma preferirei avere tutto in R in modo che la procedura sia completamente riproducibile e in modo da poter codificare dozzine di variazioni automaticamente .

Ho biforcato il codice Rcartogram su github e ho aggiunto i miei sforzi fino ad ora here.

In sostanza ciò che questa demo fa è creare uno SpatialGrid sulla mappa, cercare la densità di popolazione in ogni punto della griglia e convertirla in una matrice di densità nel formato richiesto per il cartogram() su cui lavorare. Fin qui tutto bene.

Ma, come interpolare i punti della mappa originale in base all'uscita di cartogram()?

Ci sono due problemi qui. Il primo è quello di ottenere la mappa e la griglia nelle stesse unità per consentire l'interpolazione. Il secondo è accedere a ogni punto di ogni poligono, interpolarlo e tenerli tutti nel giusto ordine.

La griglia è in unità di griglia e la mappa è in unità proiettate (nel caso dell'esempio longlat). O la griglia deve essere proiettata in longlat o la mappa in unità di griglia. Il mio pensiero è di creare un CRS falso e usarlo insieme alla funzione spTransform() in package(rgdal), poiché questo gestisce ogni punto dell'oggetto con il minimo sforzo.

L'accesso a ogni punto è difficile perché sono diversi livelli nell'oggetto SpPDF: oggetto> poligoni> poligoni> linee> coords, penso. Qualche idea su come accedervi mantenendo intatta la struttura della mappa generale?

Answer 1

Questo problema può essere risolto con il pacchetto getcartr, disponibile su Chris Brunsdon's GitHub, come splendidamente spiegato nel post del blog this.

La funzione quick.carto fa esattamente ciò che si desidera - prende uno SpatialPolygonsDataFrame come input e ha un SpatialPolygonsDataFrame come output.

Riproducendo l'essenza della esempio nel post qui nel caso in cui il link si esaurisce, con il mio stile mescolato in & errori di battitura fisse:

(Shapefile; World Bank population data)

library(getcartr) 
library(maptools) 
library(data.table) 

world <- readShapePoly("TM_WORLD_BORDERS-0.3.shp") 
#I use data.table, see blog post if you want a base approach; 
# data.table wonks may be struck by the following step as seeming odd; 
# see here: http://stackoverflow.com/questions/32380338 
# and here: https://github.com/Rdatatable/data.table/issues/1310 
# for some background on what's going on. 
world@data <- setDT(world@data) 

world.pop <- fread("sp.pop.totl_Indicator_en_csv_v2.csv", 
        select = c("Country Code", "2013"), 
        col.names = c("ISO3", "pop")) 

world@data[world.pop, Population := as.numeric(i.pop), on = "ISO3"] 

#calling quick.carto has internal calls to the 
# necessary functions from Rcartogram 
world.carto <- quick.carto(world, world$Population, blur = 0) 

#plotting with a color scale 
x <- world@data[!is.na(Population), log10(Population)] 
ramp <- colorRampPalette(c("navy", "deepskyblue"))(21L) 
xseq <- seq(from = min(x), to = max(x), length.out = 21L) 
#annoying to deal with NAs... 
cols <- ramp[sapply(x, function(y) 
    if (length(z <- which.min(abs(xseq - y)))) z else NA)] 

plot(world.carto, col = cols, 
    main = paste0("Cartogram of the World's", 
        " Population by Country (2013)"))