recogeoパッケージで市町村合併に伴う再集計をやってみる

市区町村単位で集計された統計データで時点間の比較をしようとすると、市町村合併によって2時点間で地域の対応が取れないことがあります。必要に応じて自前でスクリプトを書いて対応していますが、よさげなパッケージ(recogeo package)がリリースされていたのでちょっと試してみました。

2005年と2010年の国勢調査で得られた福岡県の人口データを例にやってみます。データは統計GISの小地域境界データのShapeファイルを使います。

パッケージのロード

library(tidyverse)
library(sf)
library(rmapshaper)

2005年データのダウンロードと展開

url2005 <- "https://www.e-stat.go.jp/gis/statmap-search/data?dlserveyId=A002005212005&code=40&coordSys=1&format=shape&downloadType=5"

tf <- tempfile()
curl::curl_download(url2005, tf)
unzip(tf, exdir = "pop2005", junkpaths = TRUE)
unlink(tf)
rm(tf)

2005年データの読み込み

st_read()シェープファイルを読み込む。海域のポリゴンを除去するために、HCODE == 8101で町丁・字等に限定。rmapshaper::ms_dissolve()で市区町村単位でディゾルブ。市区町村コードをPREFCITYを結合して作成。rmapshaper::ms_simplify()でジオメトリを簡素化。rmapshaper::ms_filter_islands()で面積の小さい離島を除去。

pop2005 <- st_read("./pop2005", options = "ENCODING=CP932",
                   stringsAsFactors = FALSE) %>% 
  filter(HCODE == 8101) %>% 
  select(PREF, CITY, CITY_NAME, JINKO, SETAI) %>% 
  ms_dissolve("CITY_NAME", sum_fields = c("JINKO", "SETAI"),
              copy_fields = c("PREF", "CITY")) %>% 
  mutate(CITY_CODE = str_c(PREF, CITY)) %>% 
  select(-PREF, -CITY) %>% 
  ms_simplify() %>% 
  ms_filter_islands(min_area = 2000000) %>% 
  rename(JINKO2005 = JINKO, SETAI2005 = SETAI) %>% 
  st_set_crs(pop2005, 4612) %>% st_transform(2444)

2010年も同様に

url2010 <- "https://www.e-stat.go.jp/gis/statmap-search/data?dlserveyId=A002005212010&code=40&coordSys=1&format=shape&downloadType=5"
tf <- tempfile()
curl::curl_download(url2010, tf)
unzip(tf, exdir = "pop2010", junkpaths = TRUE)
unlink(tf)
rm(tf)

pop2010 <- st_read("./pop2010", options = "ENCODING=CP932",
                   stringsAsFactors = FALSE) %>% 
  filter(HCODE == 8101) %>% 
  select(PREF, CITY, CITY_NAME, JINKO, SETAI) %>% 
  ms_dissolve("CITY_NAME", sum_fields = c("JINKO", "SETAI"),
              copy_fields = c("PREF", "CITY")) %>% 
  mutate(CITY_CODE = str_c(PREF, CITY)) %>% 
  select(-PREF, -CITY) %>% 
  ms_simplify() %>% 
  ms_filter_islands(min_area = 2000000) %>% 
  rename(JINKO2010 = JINKO, SETAI2010 = SETAI) %>% 
  st_set_crs(pop2005, 4612) %>% st_transform(2444)

地図の描画

市区町村数は市町村合併により、97から72に減少しています。

nrow(pop2005)
## [1] 97
nrow(pop2010)
## [1] 72

geom_sf()で地図の描画。

gridExtra::grid.arrange(
  ggplot(pop2005) +
    geom_sf(aes(fill = log(JINKO2005)), lwd = 0.2) +
    scale_fill_gradient(low="white", high="red"),
  ggplot(pop2010) +
    geom_sf(aes(fill = log(JINKO2010)), lwd = 0.2) +
    scale_fill_gradient(low="white", high="red"),
  ncol = 2)

f:id:nonki1974:20190203181801p:plain

recogeoパッケージを使う

ここからやっと本題。recogeo::reconcileGeographies()関数を使うと、2つの空間オブジェクトに含まれるポリゴン同士の関係(同じsame / 含むcontains / 重なるintersects)の判定をしてくれます。ただ、実際にはこれらの判定がそのままではうまくいかないことがあります。例えば、若松区40103)の場合、

pop2005.wakamatsu <- pop2005 %>% filter(CITY_NAME == "若松区")
pop2010.wakamatsu <- pop2010 %>% filter(CITY_NAME == "若松区")

ggplot() +
  geom_sf(data = pop2005.wakamatsu, color = "red", fill = NA) +
  geom_sf(data = pop2010.wakamatsu, color = "blue", fill = NA)

f:id:nonki1974:20190203191056p:plain

このように微妙にずれていますので、実際には同じ自治体なのですが、いずれも包含関係になりません。

そこで片方について、500mバッファをとってみます。

ggplot() +
  geom_sf(data = pop2005.wakamatsu, color = "red", fill = NA) +
  geom_sf(data = st_buffer(pop2010.wakamatsu, 500), color = "blue", fill = NA)

f:id:nonki1974:20190203191612p:plain

これなら確実に2010年のポリゴンで2005年のポリゴンを含むことができます。逆も同様に実行した場合に、包含関係が成立すれば、2つのポリゴンは同一自治体のものと判定するという感じになっている(と推測してます)。一方の包含関係のみが成立する場合はcontainになって、今回のデータの場合は、合併があったことを示すものになります。

次に、intersectについてですが、そのままintersectかどうかのチェックをすると、多くの隣接自治体がintersectと判定されてしまいます。2005年の若松区は2010年の若松区芦屋町と八幡西区とintersectの関係にあると判定されます。

st_intersects(pop2005.wakamatsu, pop2010)
## 1: 3, 8, 49

ggplot() +
  geom_sf(data = pop2010[c(3, 8, 49),], color = "blue", fill = NA) +
  geom_sf(data = pop2010.wakamatsu, color = "red", fill = NA)

f:id:nonki1974:20190203193921p:plain

これを抑止するために、intersectしたときの共通部分の面積が、一定以上の値でないと、intersectしているとみなさないということにします。今回のデータでは合併のみしか起こっていませんので、intersectの判定は生じないはずです。つまり、intersectが出現しないように、その閾値を設定すればよいことになります。

実際の判定に、recogeo::reconcileGeographies()関数を使います。

devtools::install_github('fraba/recogeo')
library(recogeo)

pop_recogeo <- reconcileGeographies(
  pop2010, pop2005, "CITY_CODE", "CITY_CODE",
  dist_buffer = 500, min_inters_area = 30000)

これで、pop_recogeoにこの期間での合併情報がdata.frameの形式で格納されます。以下は、一部の出力ですが、2006年にあった飯塚市40205)と周辺の4自治体の合併情報が確認できます。

pop_recogeo %>% arrange(unigeokey_A) %>% slice(18:24)
##   unigeokey_A unigeokey_B       type
## 1       40205       40205 AcontainsB
## 2       40205       40425 AcontainsB
## 3       40205       40426 AcontainsB
## 4       40205       40427 AcontainsB
## 5       40205       40428 AcontainsB
## 6       40206       40206       same
## 7       40207       40207       same

この情報を使って、recogeo::reconcileData()関数で人口と世帯数の集計をしてみます。

pop_recodat <- 
  reconcileData(pop_recogeo, pop2010, pop2005,
                "CITY_CODE", "CITY_CODE",
                varA = c("JINKO2010", "SETAI2010"),
                varB = c("JINKO2005", "SETAI2005"),
                return_spatial = "A")

pop_recodat %>% as_tibble() %>% select(1:5) %>% slice(1:5)
## # A tibble: 5 x 5
##   .unigeokey_new JINKO2010 SETAI2010 JINKO2005 SETAI2005
##           <dbl>     <dbl>     <dbl>     <dbl>     <dbl>
## 1              1    104469     44232    108677     44364
## 2              2     85167     33495     87340     32754
## 3              3     61583     28064     63714     28417
## 4              4    181936     89036    183286     87459
## 5              5    214793     86420    214624     83716

2010年の行政区域で2005年の人口と世帯数が再集計されています。 これで、世帯当たりの人数を地図上に可視化してみると次のようになります。

pop_recodat %>% 
  mutate(`2010` = JINKO2010 / SETAI2010,
         `2005` = JINKO2005 / SETAI2005) %>% 
  select(`2005`, `2010`) %>%
  gather(year, n_household, `2005`, `2010`) %>%
  ggplot(aes(fill = n_household)) +
  geom_sf() +
  facet_wrap(.~year) +
  scale_fill_gradient(low="white", high="red")

f:id:nonki1974:20190203203028p:plain

とりあえず、合併のみという単純な場合なら使えそうな感じです。政令指定都市への移行などで生じる分割や、もっと複雑な場合もありますが、その検証はまた別の機会に。

Enjoy!

TinyTeXのインストール & RmarkdownでPDF on Ubuntu 18.04 LTS

TinyTex のインストール

Motivation

TeXのインストールはヘビーなので、なるべく簡素化したい。TeX Liveとか数GBをダウンロードせんといかん。

環境

  • Ubuntu 18.04 LTS on Hyper-V
  • 言語パックインストール済み

Atomのインストール

TeXのエディタとしてAtomを使う。

sudo add-apt-repository -n -y ppa:webupd8team/atom
sudo apt update
sudo apt install atom

パッケージのインストール

apm install japanese-menu latex language-latex latexer pdf-view

TinyTexのインストール

wget -qO- "https://yihui.name/gh/tinytex/tools/install-unx.sh" | sh

TinyTeXのインストール(for R users)

すでにRが入っていれば、以下でOK。

install.packages('tinytex')
tinytex::install_tinytex()

PATHの設定

export PATH=$HOME/.TinyTeX/bin/x86_64-linux:$PATH

.bashrcに追加しとく

TeXの日本語環境をインストール

# リポジトリの設定
tlmgr option repository http://ftp.jaist.ac.jp/pub/CTAN/systems/texlive/tlnet/
# collection-langjapaneseのインストール
tlmgr install collection-langjapanese

Notoフォントの設定

Ubuntuの日本語言語パックでNotoフォントが/usr/share/fonts/opentype/notoに入っているので、それを使う。入っていない場合はインストールする。

mkdir texmf-local
cd texmf-local
mkdir fonts
cd fonts
ln -s /usr/share/fonts/opentype/ opentype
mktexlsr

Atomの設定とタイプセットのテスト

Atomlatexパッケージの設定でEngineuplatexにする。

\documentclass[a4j, uplatex]{jsarticle}

\usepackage[noto-otc]{pxchfon}

\begin{document}
  \section{はじめに}
  $X \sim N(\mu, \sigma^2)${\bf 確率密度関数}\[
  f(x) = \frac{1}{\sqrt{2\pi}\sigma}
  \exp\left\{-\frac{(x-\mu)^2}{2\sigma^2}\right\}
  \]
  である。
\end{document}

完成!

f:id:nonki1974:20181227002409p:plain
AtomTeX

knitrでPDF出力

この流れでRmarkdownからpdf-documentとかbeamer presentationの作成までやってしまいたい。

Rのインストール

/etc/apt/sources.listに以下を追加

deb https://cloud.r-project.org/bin/linux/ubuntu bionic-cran35/

キーの登録、パッケージリスト更新、インストール

sudo apt-key adv --keyserver keyserver.ubuntu.com --recv-keys E298A3A825C0D65DFD57CBB651716619E084DAB9
sudo apt update
sudo apt install r-base r-base-dev

RStudioのインストール

wget https://download1.rstudio.org/rstudio-xenial-1.1.463-amd64.deb
sudo apt install gdebi-core
sudo gdebi rstudio-xenial-1.1.463-amd64.deb

IPAexフォント入れとく

sudo apt install fonts-ipaexfont

Rmdファイルの作成

---
title: "TinyTeXを使ってみよう"
author: "nonki1974"
date: "2018年12月26日"
output:
  pdf_document:
    latex_engine: xelatex
documentclass: bxjsarticle
classoption: xelatex, ja=standard
geometry: no
header-includes:
  - \usepackage{zxjatype}
  - \setjamainfont{Noto Serif CJK JP}
  - \setjamonofont{Noto Sans Mono CJK JP}
  - \setjasansfont{Noto Sans CJK JP Medium}
  - \setmainfont{Noto Serif CJK JP}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
knitr::opts_chunk$set(dev="cairo_pdf",
                      dev.args=list(family="Noto Sans CJK JP"))
```

## R Markdownのテスト
`pressure`データのプロットです。
```{r}
plot(pressure, xlab = "温度[℃]", ylab = "蒸気圧[mm]")
```

knitする

できた〜

f:id:nonki1974:20181227010200p:plain
knitrでPDF

宿題

Windowsでも動作確認

sfパッケージ&国土数値情報でコロプレス図

Rで地図を描画する際はspパッケージが必須だったわけですが、今後sfパッケージの方が便利にになるらしいです。dplyrと相性がいいんでしょうね。このでやっている、ごみ排出量のコロプレス図の作成をsfパッケージでやってみることにします。

シェープファイルの取得

国土数値情報からシェープファイルをダウンロード。kokudocuuchiパッケージを使って、URLを確認。

library(dplyr)
kokudosuuchi::getKSJSummary() %>% filter(title=="行政区域")
## # A tibble: 1 x 5
##   identifier    title   field1 field2 areaType
##        <chr>    <chr>    <chr>  <chr>    <chr>
## 1        N03 行政区域 政策区域      -        3

ということなので、福岡県(prefCode = 40)を指定してURL一覧を取得。そこから2015年のもののみ取り出す。

kokudosuuchi::getKSJURL("N03", prefCode = 40) %>%
  filter(year==2015) %>%  select(zipFileUrl) %>%
  pull(1) -> url

ファイルをダウンロードしてカレントディレクトリに展開。

tf <- tempfile()
curl::curl_download(url, tf)
unzip(tf, exdir = getwd(), junkpaths = TRUE)
unlink(tf)
rm(tf)

シェープファイルの読み込み→加工→白地図の描画

sfパッケージをロードして、シェープファイルを読み込む。

library(sf)
d <- st_read(".", options = "ENCODING=CP932", stringsAsFactors = FALSE)

政令指定都市については、ごみ排出量のデータは、全体での集計になっているのに対して、シェープファイルは区ごとのデータになっているため、政令指定都市全体をaggregate()関数を使って1つのポリゴンにまとめる。まとめた後に、市区町村名(N03_004)と、市区町村コード(N03_007)の修正が必要。

d %>% filter(N03_003 %in% c("北九州市","福岡市")) %>%
  aggregate(by = list(.$N03_003), FUN = "head", n=1) %>%
  select(-1) %>% 
  mutate(N03_004 = N03_003,
         N03_007 = case_when(N03_004 == "北九州市" ~ "40100",
                             N03_004 == "福岡市" ~ "40130")) -> d.seirei

他の市区町村についても、国土数値情報のデータは1つの市町村が多くのポリゴンの集まりになっているので、これらも1つのポリゴンにまとめる。

d %>% filter(!(N03_003 %in% c("北九州市", "福岡市"))) %>%
  aggregate(by = list(.$N03_007), FUN = "head", n=1) %>%
  select(-1) -> d.others

これらをまとめる。

rbind(d.seirei, d.others) %>% arrange(N03_007) -> d.new

白地図を出力。

plot(st_geometry(d.new), border="grey", col="white")

f:id:nonki1974:20170712081741j:plain

ごみ排出量データの取得と、sfオブジェクトへのマージ

ダウンロードするファイルのURLと保存先ファイル名の指定。e-statのAPI使ってやると便利かもしれないけど、idの取得が必要になるので今回はパス。

waste_url <- "http://www.env.go.jp/recycle/waste_tech/ippan/h27/data/shori/city/40/01.xls"
waste_fn <- "waste_fuk_h27.xls"

ダウンロードして、readxlパッケージでExcelのファイルを読み込む。列名を設定して、特に使う項目は分かりやすい列名に変更。

curl::curl_download(waste_url, waste_fn)
waste_fuk <- readxl::read_excel(waste_fn, sheet = 1, skip = 6)
colnames(waste_fuk) <- paste0("V", seq(1:ncol(waste_fuk)))
waste_fuk %>% rename(jcode=V2, name=V3, waste_total=V12) -> waste_fuk

sfオブジェクトへマージ。

left_join(d.new, waste_fuk, by = c("N03_007" = "jcode")) -> d.new2

コロプレス図(塗り分け地図)の作成。

library(classInt)
library(RColorBrewer)
ci <- classIntervals(d.new2$waste_total, style = "pretty")
d.new2 %>% select(waste_total) %>%
  plot(col=findColours(ci, rev(brewer.pal(7, "RdYlGn"))))

f:id:nonki1974:20170712082011j:plain

ggplot2のgeom_sf()を使った地図の描画

開発版のggplot2では、geom_sf()を使って、sfオブジェクトからggplot2形式のプロットが作成できます。 Windowsの場合は開発ツールをインストールしてから、以下のようにすればインストールできます。

devtools::dev_mode()
devtools::install_github("hadley/scales")
devtools::install_github("tidyverse/ggplot2")
library(ggplot2)
ggplot(d.new2, aes(fill=waste_total)) + geom_sf() +
  scale_fill_distiller(palette="RdYlGn") + theme_bw()

f:id:nonki1974:20170712082145j:plain

enjoy!

vdmR package(Tokyo.R#55のLT)

Tokyo.R#55に参加してきました。LTで自分の作成したvdmR packageを紹介しましたが、スライドのほとんどがスクリーンキャストだったので、スライドではなくこちらにまとめておきます。

vdmR packageとは

  • VDM:Visual Data Mining
  • 探索的データ解析のための対話的統計グラフィックス
  • 最新バージョンは0.2.2(最初のリリースは2014-09-20)

基本機能

  • データフレームから複数のプロットを出力し、1つのプロットにおけるデータ点の選択によって、他のプロットの対応する点やデータテーブルの対応するレコードがハイライトする(Multiple Linked View
  • プロットはJavaScriptを含んだHTMLやSVGで出力され、ブラウザ上に表示される。

デモ1(基本機能)

デモ2(コロプレス図)

デモ1のコード

# パッケージのインストールとロード
install.packages("vdmR")
library(vdmR)

# サンプルデータの読み込み
data(vsfuk2012)
head(vsfuk2012[,1:5])

# 死亡率と出生率の散布図
# Working DirectoryにHTMLとSVGが出力される
vscat(MortalityRate, FertilityRate, vsfuk2012,
      "scat01", "vsfuk2012", color=Type, size=pop_male)

# 婚姻率のヒストグラム
# Working DirectoryにHTMLとSVGが出力される
vhist(MarriageRate, vsfuk2012, "hist01", "vsfuk2012",
      fill=I("darkgreen"), color=I("black"))

# Working DirectoryにデータテーブルのHTMLが出力され
# デフォルトブラウザで表示される
vlaunch(vsfuk2012, "main", "vsfuk2012")

# inline表示
vlaunch(vsfuk2012, "main", "vsfuk2012", iframe=T)

デモ2のコード

# 平行座標プロットの出力
vpcp(vsfuk2012, 4:17, "pcp1", "vsfuk2012",
     groupColumn="Type", scale="uniminmax", missing="min10")

# コロプレス図

# maptoolsの読み込み(事前にインストール)
library(maptools)

# サンプルで同梱されているシェープファイルのパスを指定
shp.path <- file.path(system.file(package="vdmR"),
                      "etc/shapes/fukuoka2012.shp")

# シェープファイルの読み込み
vsfuk2012.spdf <- readShapeSpatial(shp.path, IDvar="CityCode")

# 属性テーブルのプレビュー
head(vsfuk2012.spdf@data)

# 塗り分けのための色の設定
frcol <- ggplot2::scale_fill_gradient2(low="blue", mid="white", high="red",
                                       midpoint=median(vsfuk2012$FertilityRate))

# コロプレス図の出力
vcmap(shp.path, vsfuk2012, "CityCode", "CityCode", "map1", "vsfuk2012",
      fill=FertilityRate, ggscale=frcol)

# データテーブルとプロットの表示(インライン)
vlaunch(vsfuk2012, "main", "vsfuk2012", iframe=T)

2014年のパッケージダウンロード数ランキング(CRAN rstudio mirror)

ちょっと調べる必要があったので,ついでにメモ.

installrパッケージで,CRAN rstudio mirrorのダウンロードログを取得して集計(上位20パッケージ).まあ,rstudio mirrorなので,全体からすると多少偏りはあると思います.

f:id:nonki1974:20151120182152p:plain

        package download
1       ggplot2   880791
2          Rcpp   815665
3        BayHaz   771257
4          plyr   767302
5       stringr   677974
6        digest   676556
7          clhs   663409
8  RColorBrewer   609341
9      reshape2   605751
10   colorspace   540774
11       scales   508104
12        rJava   502119
13     labeling   476075
14        proto   463113
15      munsell   457374
16       gtable   447665
17    dichromat   447135
18          DBI   445992
19       bitops   427203
20          car   403301

leafletRでシェープファイルからOSM上に塗り分け地図

Leaflet - a JavaScript library for mobile-friendly maps

というインタラクティブ地図のためのJavaScriptライブラリがある。GeoJSONファイルを用意して、leafletのスクリプトを書いてやれば、OpenStreetMapの上にインタラクティブ機能付きでデータを可視化できる。

これの作業をR上でやってくれるのがleafletRパッケージ。

福岡県の市町村別高齢化率(平成26年)の塗り分け地図(コロプレス図)を作成してみる。

# leafletRパッケージを使って、Shapefileから
# leafletJSを使ったインタラクティブな地図(コロプレス図)を作成
# leafletJS: http://leafletjs.com

library(leafletR)
library(rgdal)

tmp <- tempdir()

# 国土数値情報の行政界(平成26年度、福岡県)のシェープファイルを編集して、
# 各市町村の高齢化率を属性テーブルに追加した。
# 高齢化率は福岡県の高齢者人口及び高齢化率の推移
# (http://www.pref.fukuoka.lg.jp/contents/koureika-suii26-2.html)
# のエクセルファイルより取得。
url <- "http://stat.fwu.ac.jp/~fujino/hateblo/20150510/fukuoka.zip"

file <- basename(url)
download.file(url, file)
unzip(file, exdir=tmp)

# シェープファイルをSpatialPolygonsDataFrameとして取り込む
fukuoka <- readOGR(dsn=paste0(tmp,"/fukuoka"), 
                        layer="fukuoka", encoding="Shift_JIS")

# GeoJSONオブジェクトに変換
# ちょっと時間がかかる
fukuoka.dat <- toGeoJSON(data=fukuoka, name="fukuoka", dest=tmp)

# スタイルの設定
# propは塗り分ける属性テーブルの列名、breaksは階級、
# style.valは階級の色、legは凡例のタイトル
fukuoka.style <- styleGrad(prop="agingrate", breaks=seq(10,40,by=5),
                    style.val=rev(heat.colors(6)), leg="Aging Rate",
                    fill.alpha=0.7, lwd=0.5)

# leafletJSのスクリプトを含むHTMLファイルを生成
# シェープファイルの属性テーブルにマルチバイト文字が含まれていて
# 文字コードがUTF-8でない場合にはここでエラーが出る
fukuoka.map <- leaflet(data=fukuoka.dat, dest=tmp, title="Aging Rate in Fukuoka",
                    base.map="osm", style=fukuoka.style, popup="*")
# ブラウザに表示
fukuoka.map

ポリゴンをクリックすると属性テーブルの情報がポップアップで表示される。
f:id:nonki1974:20150510225238p:plain

動くのは以下のリンクで。
Aging Rate in Fukuoka

githubで公開されているRStudioのleafletパッケージの方がポップアップのカスタマイズができたりと高機能な感じはするけど、まだ試していない。

gridSVGで散布図上の点のIdentification (2)

前回のエントリ

のサンプルを確認してみたら最新の環境ではうまくいかなかったので、アップデート。

grid.garnish("geom_point.points.1", onclick="info(evt)")

で点をクリックしたときのイベントを登録しても、出力のSVGに全く反映されない。とりあえずは以下のようにすれば、期待する結果が得られる。

grid.set("geom_point.points",
         garnishGrob(grid.get("geom_point.points"), onclick="info(evt)"))

grid.garnish()のオプションでも何とかなるかもしれないが、いまのところ原因は不明。ggplot2が出力するglobの構造が変わってる可能性があるかも。

もう1つはJavaScript

id = e.target.correspondingUseElement.id;

chromeの最新版でcorrespondingUseElementがundefinedになってしまう。

id = e.target.id;

でうまくいく。単に使わなくすればよいだけ。以前のエントリでchromeのみで動作と書いたが、結果的にこれでfirefoxでもOKになった!でも当面はクロスブラウザ対策が必要かな。

全体のコードは以下。

pdf(file=NULL, width=7, height=5)

ggplot(iris, aes(Petal.Width, Petal.Length)) + 
  geom_point(aes(colour=Species, size=Sepal.Width))

grid.script(paste("var x=", toJSON(iris$Petal.Width),";",sep=""))
grid.script(paste("var y=", toJSON(iris$Petal.Length),";",sep=""))

grid.script("
	info = function(e){
		id = e.target.id;
		id = parseInt(id.replace('geom_point.points.1.',''));
		alert('Petal.Width='+x[id-1]+'\\nPetal.Length='+y[id-1]);
	}
")

grid.gedit("geom_point.points", name="geom_point.points")
grid.set("geom_point.points",
         garnishGrob(grid.get("geom_point.points"), onclick="info(evt)"))

grid.export("iris_plot.svg")
dev.off()