7

私はこの質問で尋ねられたのと同じことをしようとしています. Cartogram + choropleth map in Rですが、 SpatialPolygonsDataFrame から始めて、同じタイプのオブジェクトになることを望んでいます。

オブジェクトをシェープファイルとして保存し、scapetoadを使用して再度開き、元に戻すこともできますが、手順を完全に再現できるように、また数十のバリエーションを自動的にコーディングできるように、すべてを R 内に置きたいと思います。

私は github で Rcartogram コードをフォークし、これまでの努力をここに追加しました。

基本的に、このデモが行うことは、マップ上に SpatialGrid を作成し、グリッドの各ポイントで人口密度を検索し、これをcartogram()が作業するために必要な形式の密度マトリックスに変換することです。ここまでは順調ですね。

しかし、の出力に基づいて元のマップ ポイントを補間する方法はcartogram()?

ここには 2 つの問題があります。1 つ目は、マップとグリッドを同じ単位にして、補間を可能にすることです。2 つ目は、すべてのポリゴンのすべてのポイントにアクセスして補間し、それらをすべて正しい順序に保つことです。

グリッドはグリッド単位で、マップは投影単位です (例の longlat の場合)。グリッドを longlat に投影するか、マップをグリッド単位に投影する必要があります。spTransform()私の考えでは、偽の CRS を作成し、これを の関数と一緒に使用します。これpackage(rgdal)は、最小限の手間でオブジェクト内のすべてのポイントを処理するためです。

SpPDFオブジェクトにいくつかのレイヤーがあるため、すべてのポイントにアクセスすることは困難です。オブジェクト>ポリゴン>ポリゴン>ライン>座標だと思います。マップ全体の構造をそのまま維持しながら、これらにアクセスする方法はありますか?

4

1 に答える 1

2

この問題は、このブログ投稿で美しく説明されているように、 Chris Brunsdon の GitHubgetcartrで入手できるパッケージで解決できます。

quick.carto関数はまさにあなたが望むことを行います - 入力として a を取り、出力としてSpatialPolygonsDataFramea を持ちSpatialPolygonsDataFrameます。

リンクが機能しなくなった場合に備えて、ブログ投稿の例の本質を再現し、私自身のスタイルを混ぜて誤植を修正します。

(シェープファイル;世界銀行の人口データ)

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)"))

ここに画像の説明を入力

于 2015-09-08T23:43:52.260 に答える