文系のための「地理情報標準」（１）

様々な情報管理システムの中でも、特に空間情報を扱うシステムのことを、

地理情報システム（GIS: Geographic Information Syste）あるいは、

地理空間情報システム（GIS: Geospatial Information System）と呼ぶ。

略称を使うとどちらも同じで、略称でGISと呼ばれているものがある。

教科書的な定義では、「地理空間情報を取得、保存、統合、管理、分析、

伝達して、空間的意思決定を支援するコンピュータベースの技術」である。

（村山裕司・柴崎亮介（騙）(2008). 『シリーズGIS第一巻　GISの理論』, 朝倉書店.）

おそらく、このGISというのが一体どのようなシステムなのか、

今となっては、当たり前すぎる技術となっていて、

明確な定義をすること自体が困難なものとなっている。

良い言い方ではないけれど「高度な電子地図システム」としておく。

実を言うと、このGISの歴史を述べるのは少々厄介で、

学界、政界、財界、のそれぞれの視点という思惑によって、

その起源は異なったコンテキストで語られている。

今では、GISのたとえ話で頻繁に登場するカーナビやGoogle Mapなど、

かつては、GISとは異なるものとして扱われていたし、

一緒にすることを極端に避けようする人も少なかった…。

いやいや、語り始めると本題に入れないので、今は、この話は置いておく。

そうそう、ここで議論するべきことは、GISの歴史や定義の話ではなく、

もっと技術的な内容である。すなわち、GISのデータ管理の話で、

特に、その基礎的な設計方法について整理するのであった。

GISにおけるデータを設計する上で一番厄介な問題は、

GISが現実世界のありのままの姿に対する人の認識、

というのが基礎的な概念として必要となっている点である。

つまり、「哲学的な話」と「技術的な話」の両方が必要となってくる。

実際に、GISのデータに関わる国際標準『地理情報標準』は、

かなり抽象的な認識論に基づいていて、多くの技術者を悩ましている。

おそらく、この国際標準を理解できる素養があるのは、

高度な技術力を持つ理工系エンジニアではなく、高い教養と

観念論的な思考が可能な文系出身のエンジニアであると思う。

というのは、GISでは対象となる現実世界をコンピュータ上に

モデル化する必要があって、したがって、現実世界の捉え方

というのが、どうしても避けられない問題となってくる。

そうした捉え方の一つが「オブジェクト指向GIS（OOGIS）」であり、
『地理情報標準』という標準は、このOOGISの考え方に基づいているのだが、
まずは、この標準について、簡単に整理しておくことにする。

そもそも、先ほどから『地理情報標準』と言っているけれど、
国際的には「ISO 19100」シリーズと呼ばれていて、
国際標準化機構という機関が国際標準として策定している。

この国際標準化機構は、ISOと言う略称で知られていて、

生産品、製品、素材、材料、プロセスに関する規格が中心に扱われる。
有名なものに、クレジットカードの寸法やネジの規格などがある。

ISOでは、技術委員会の下で標準化が進められていて、
ISO 19100シリーズの場合は、211番目の技術委員会「ISO/TC 211」で
検討が重ねられている。

この委員会での議論の範囲は、「地球上の場所に直接的、
あるいは間接的に関連付けられた事象の情報について、
体系的な一連の標準として確立すること」である。

また、具体的には、「異なるユーザーやシステム間で電子化された
データを管理、取得、加工、分析、接続、表現、伝達するための
方法論、ツール、サービスを明確化する」ことを目的としている。

何やら、難しい話に聞こえるが、要するに、地理的な情報に
関わる情報の取扱全般におよぶプロセスを整理するための
国際標準を作るのがISO/TC 211の目的となっている。

ISO/TC 211では、さらに、議論するべき内容によって、
ワーキング・グループが置かれていて、
それぞれのワーキング・グループで標準化を行っている。

• WG 1: フレームワークと参照モデル
• ISO 19101: 参照モデル
• ISO 19103: 概念スキーマ言語
• ISO 19104: 専門用語
• ISO 19105:  適合性と検証
• ISO 19121: 画像とグリッドデータ
• WG 2: 地理空間データモデルと操作（解散）
• ISO 19107: 空間スキーマ
• ISO 19108: 時間スキーマ
• ISO 19109: アプリケーション・スキーマの規則
• ISO 19123: 被覆の幾何と機能のためのスキーマ
• WG 3: 地理空間データの管理（解散）
• ISO 19110: 地物カタログの方法論
• ISO 19112: 地理識別子による空間参照
• WG 4: 地理空間サービス
• ISO 19116: 位置情報サービス
• ISO 19117: 描画法
• ISO 19118: 符号化
• ISO 19119: サービス
• ISO 19125: シンプル・フィーチャー・アクセス
• ISO 19128: Webマップ・サーバーのインタフェース
• ISO 19136: Geography Markup Language（GML）
• ISO 19142: Webフィーチャー・アクセス
• ISO 19149: 地理情報のための権利表示言語（GeoREL）
• ISO 19161: 測地系
• ISO 19164: 登録サービス
• WG 5: プロファイルと機能に関する標準（解散）
• ISO 19106: プロファイル
• WG 6: 画像
• ISO 19120: 機能的標準
• ISO 19129: 画像、グリッド、被覆の枠組み
• ISO 19130: 地理情報取得のための画像センサーのモデル
• ISO 19159: リモートセンシング・センサーのキャリブレーション
• ISO 19163: 画像とグリッドデータの内容コンポーネントと符号化規則
• WG 7: 情報の公共性
• ISO 19115: メタデータ
• ISO 19122: 人材の資格と認定
• ISO 19126: フィーチャーの概念辞書と登録
• ISO 19137: 空間スキーマのためのコア・プロファイル
• ISO 19139: メタデータの実装
• ISO 19144: 分類システム
• ISO 19150: オントロジー
• ISO 19152: 土地管理ドメイン（LADM）
• ISO 19160: 住所表記
• ISO 19165: デジタルデータとメタデータの保存
• WG 8: ロケーション・ベース・サービス（解散）
• ISO 19132: ロケーション・ベース・サービスの参照モデル
• ISO 19133: ロケーション・ベース・サービスの追跡とナビゲーション
• ISO 19134: ロケーション・ベース・サービスのマルチモーダルルーティング
• WG 9: 情報管理
• ISO 19111: 座標による空間参照
• ISO 19127: 測地系コードとパラメータ
• ISO 19131: データの製品仕様書
• ISO 19135: 項目登録のための手順
• ISO 19145: 地点の表現の登録
• ISO 19146: クロス・ドメインの語彙
• ISO 19153: 地理空間の電子的な権利管理参照モデル（GeoDRM RM）
• ISO 19156: 観察と計測
• ISO 19157: データ品質
• ISO 19158: データ提供の品質の確実性
• ISO 19162: WKTによる座標参照系の表現
• WG 10: ユビキタスと公共アクセス
• ISO 19147: 移送ノード
• ISO 19148: 線形参照
• ISO 19154: ユビキタス公共アクセス参照モデル

実際には、同じ標準を異なるワーキング・グループで議論されていたり、

審議内容が別のワーキング・グループに移管されたり、

あるいは、解散したワーキング・グループも存在するので、

ここに書かれている通りでは無い部分もある。

そういった、細かい点は置いておいても、これを見てみると解るように、

様々な内容について議論されていて、各ワーキング・グループで、

関連する標準について議論されていることがわかる。

残念ながら、私自身、全ての標準の中身を理解できている訳ではないし、

これらの標準の全てを理解できている人は、世界でも居るのか…。

ところで、本ブログでは、ISO 19100 シリーズとは呼ばずに、

『地理情報標準』という言葉を使っているが、

その理由に関しても、少し、説明しておく必要がある。

実は、国際的な標準としてはISO 19100として通用しているけれど、

日本国内に同標準を適用するには、日本にのみ関わる部分や、

日本の特有の事情で補足あるいは追加するべき点がある。

そういった部分を精査して、国内向けに別の標準が作られている。

その標準は、日本工業規格（JIS）のJIS X7100シリーズとして、

国内での普及が始まっている。

さらに、国土地理院は、JISX7100シリーズのうち、

標準の中でも実用上必要な部分を取り出して体系化し、

『地理情報標準プロファイル（JPGIS）』を策定している。

これとは別に、かつては、国土地理院と国内民間企業が共同研究として、

ISO/TC 211の議論に準拠して国内への適合を検討していた時期があって、

その成果が『地理情報標準（JSGI）』としてまとめられている。

細かい内容を見てみると、当然のことながら、それぞれに相違はあるけれど、

全体的な体系や考え方は同じなので、総称して『地理情報標準』としている。

この呼び方の方が、内容がよく解るというメリットもある。

とりあえず、『地理情報標準』については、この辺りで止めておいて、

少しずつ、同標準の考え方というか、OOGISの考え方について、

次回から、理解を進めていきたい。

自分のための「Rで重心ボロノイ基準？」

先週、谷岡氏と二人で相談して毎週「意味不明」な分析手法を開発することにした。
今週の課題は、「重心ボロノイ」を用いた方法の検討。
この手法は、ボロノイセルのサンプルを、そのセルの重心へと移動し続けるというもの。

逐次計算を行うので、どこで止めるかが問題になるが、
現段階では、重心の移動量がある閾値よりも小さくなったときに終了。
距離の計算方法も色々と設定できそうではあるが、今は考慮しない。

この方法、何に役立つのか？

ふむ。最終的に得られる重心ボロノイは、同じ条件下での領域だから、
つまりは、地理空間に適用した場合、最終的に得られた重心ボロノイは、
潜在的に存在する地理的要因を取り除いたものになるのか？？？

例えば、土地起伏、植生、河川などの水場、「不吉」な空間...etc
そういった地理的要因を取り除いた空間におけるボロノイってことか？
いや、よく解らない。まぁ、これで論文書くわけじゃないし...まぁ、いいや。

実用化するには、もう少し、色々と検討する余地があるが、
重心ボロノイの面積（の逆数）や、初期ボロノイのサンプル点からの移動量、
そういった値の分散を取るのも面白い。潜在的な要因を探ることができそう。

今度、何かの研究で使ってみるかな？

今回は、Rでコロプレスを描く方法も利用する。
長いソースコードをコピペすると、時々、エラーが出るのでスクリプトファイルに。
とりあえず、スクリプトファイルは以下の場所からダウンロード可能。

https://docs.google.com/open?id=0B9OtIs5BjRVhQXRaUk5ZOUEyUXM

とにかく、実行コード

---

# 最初にR スクリプトのソースを読み込む。
source("path/to/source/cis.geom.r")

# シェープファイルを読み込む
shp <- readShapePoints("path/to/shapefile/shapes.shp")

# 重心ボロノイを描く。
cv <- cis.centroidalVoronoi(shp)

# 結果をプロットする。
par(mfrow=c(1,2))

plot(deldir(shp@coords[,1], shp@coords[,2]),wlines="tess", pch=16, main="Original")

cis.choropleth(cv, "a", main="Centroidal Voronoi")

っで、Rスクリプに含まれているソースコード。

---

cis.centroidalVoronoi <- function(shp, eps=0.05, verbose=0){
# Load required packages.
require(voronoi)
require(maptools)
require(utils)

# Build original voronoi diagram.
obj <- deldir(shp@coords[,1], shp@coords[,2], verbose=0)

# Build centroidal voronoi.
c.obj <- centroidal(obj, eps=eps)

n <- length(tile.list(c.obj$T))
l <- list()
a <- vector(length=n)

# Create a text progressbar.
pb <- txtProgressBar(min=0, max=n)
for(i in 1:n){
 # Get coordinates of the tesseract.
 x <- c(tile.list(c.obj$T)[[i]]$x, tile.list(c.obj$T)[[i]]$x[1])
 y <- c(tile.list(c.obj$T)[[i]]$y, tile.list(c.obj$T)[[i]]$y[1])

 # Create a polygon from the set of coordinates.
 p <- Polygon(matrix(c(x,y),ncol=2), hole=FALSE)

 # Convert a polygon to the polygons.
 l[i] <- Polygons(list(p), i)
 a[i] <- l[[i]]@area^(-1)

 # Update the text progressbar.
 setTxtProgressBar(pb, i)
}

# Create a spatial polygon from the list of polygons.
sp <- SpatialPolygons(l)
sp <- addAttrToGeom(sp, data.frame(a), FALSE)

# Return centroidal voronoi with spatial polygons dataframe.
return(sp)
}

自分のための「Rで地理空間3D樹状図」

このコードは、シェープファイルを読み込んで、
平面上の距離から階層的クラスタリングを行い、
さらに、その結果を地理空間上に投影するというもの。

このソースコードは、私の「思い付き」に無理を言って後輩を付き合わせ、
私は、ソースコードをチェックして書き直しただけ。
同志社大学大学院文化情報学研究科の谷岡さん、長々とありがとう。

本人に許可をもらったので、このままUP。
もちろん、私も谷岡さんも一切の責任を取りません。
距離法のところは注意。Ward Method のときは、ユークリッド「平方」距離。

上記注意事項の意味が解らない人は、使用を避けた方が良いかもしれない。
他にも、細々とした所に問題があるのは理解している。
比較的、可読性は高いと思うので、見れば解るだろう。

なお、二人で相談した結果、論文クオリティに足りないということで、
このソースコードは「使い道」を知っている人に使って頂きたい。
使う時には、このブログを参考文献に...。

ソースコードを自分の用途に書き換えれば、多少は使えるかもしれない。

とりあえず、今回のサンプルデータは以下の場所に。
http://docs.google.com/open?id=0B9OtIs5BjRVhdWJmWU43V0Z3eVE

実行コードは以下の通り。

---

# Read shape file.
hoge <- readShapePoints("Krugans.shp")

# Number cluster. k <- 3
clus.method <- "average"

# Build 3D dendrogram. dendro.3d(hoge, k, clus.method)

ソースコードは以下の通り。

---

dendro.3d <- function(shp, k, c.method){
    # Read packages for this function.
    require(rgl)
    require(maptools)

    # Get coordinates of points data.
    coords <- shp@coords
    oldnames <- rownames(coords)
    rownames(coords) <- 1:nrow(coords)
    
    # Add nodes corresponding clusters.
    rowcount <- nrow(coords)

    # Conduct cluster analysis
    coords.dist <- dist(coords)
    coords.clus <- hclust(coords.dist, method=c.method)

    # Get a list of name of merge.
    coords.merge <- coords.clus$merge
    coords.merge2 <- coords.clus$merge
    coords.height <- coords.clus$height

    # Add a new column to store Z values.
    coords <- cbind(coords, rep(0, nrow(coords)))
    colnames(coords) <- c("x","y","z")

    # Change merged nodes name
    for(i in 1:nrow(coords.merge)){
            if(as.numeric(coords.merge[i,1]) > 0){
                rownum <- as.numeric(coords.merge[i,1])
                node_org <- coords.merge[rownum,1]
                node_dst <- coords.merge[rownum,2]
                node_new <- paste(node_org, node_dst, sep="")
                coords.merge[i,1] <- node_new
            }
            if(as.numeric(coords.merge[i,2] > 0)){
                rownum <- as.numeric(coords.merge[i,2])
                print(rownum)
                print(i)
                node_org <- coords.merge[rownum,1]
                node_dst <- coords.merge[rownum,2]
                node_new <- paste(node_org, node_dst, sep="")
                coords.merge[i,2] <- node_new
            }
    }

    # Combine a list of merged nodes and a list of height.
    coords.merge_and_height <- cbind(coords.merge, coords.height)

    # Rename rownames.
    rownames(coords) <- -1*as.numeric(rownames(coords))
    
    for(i in 1:(rowcount-1)){
       # print(i)
        # Get z value from height of cluster.
        z <- coords.merge_and_height[i,3]

        # Create a new node name.
        node_org <- coords.merge_and_height[i, 1]
        node_dst <- coords.merge_and_height[i, 2]
        newname <- paste(node_org, node_dst, sep="")

        # Add new column with new node name.
        coords <- rbind(coords, as.numeric(c(0,0,z)))
        rownames(coords)[i+rowcount] <- newname

        x_org <- as.numeric(coords[rownames(coords)==node_org,1])
        x_dst <- as.numeric(coords[rownames(coords)==node_dst,1])
        x_new <- (x_org + x_dst)/2

        y_org <- as.numeric(coords[rownames(coords)==node_org,2])
        y_dst <- as.numeric(coords[rownames(coords)==node_dst,2])
        y_new <- (y_org + y_dst)/2

        coords[i+rowcount, 1] <- x_new
        coords[i+rowcount, 2] <- y_new
        coords[i+rowcount, 3] <- as.numeric(z)
    }
    # Define colours for nodes.
    nodes.col <- as.numeric(cutree(coords.clus,k))+1
    nodes.marged.col <- c(nodes.col, rep(1,length(coords.height)))
    
    # Plot nodes.
    plot3d(coords,col=nodes.marged.col)
    
    # Plot branches.
    for(i in 1:(rowcount-1)){
        hoge1 <- coords[(i+nrow(shp)),]
        segments3d(rbind(coords[rownames(coords)==coords.merge[i,1],],c(coords[rownames(coords)==coords.merge[i,1],1:2],hoge1[3])),col="orange")
        segments3d(rbind(coords[rownames(coords)==coords.merge[i,2],],c(coords[rownames(coords)==coords.merge[i,2],1:2],hoge1[3])),col="orange")
        segments3d(rbind(c(coords[rownames(coords)==coords.merge[i,1],1:2],hoge1[3]),c(coords[rownames(coords)==coords.merge[i,2],1:2],hoge1[3])),col="orange")
    }
}

自分のための「RでShape file の表示」

Rでは、plot()関数やspplot()関数など、空間データの表示方法がいくつかある。
ただ、中々、痒い所に手が届かず、色々と不満。
そんな訳で、自分用の描画デバイスを作ってみようかと...。

今回のサンプルデータは以下の場所に。
https://docs.google.com/folder/d/0B9OtIs5BjRVhWm8zV1hCUFh6YnM/edit

とりあえず、表示は上手く行く。あとは、凝った表現。
どうやら、図の重ね合わせ等もできそう。
ためしに、どのようなものか挑戦してみた。やはり、若干面倒か？

以下が実行コード。

---

# Load a required library.
library(maptools)

# Read a shape file.
shp <- readShapePoly('/home/yufujimoto/Desktop/japan.shp')

# Draw polygons.
drawShape.poly(shp)

そして、以下がソースコード

---

drawShape.poly <- function(shp){
    require(grid)
 
    # Get range from bounding box.
    range.x <- c(shp@bbox[1,1],shp@bbox[1,2])
    range.y <- c(shp@bbox[2,1],shp@bbox[2,2])

    # Get widths and heights from data.
    width <- diff(range.x)
    height <- diff(range.y)

    pushViewport(
        viewport(name="layer",
                layout=grid.layout(1,1,
                       widths=width,
                       heights=height,
                       respect=TRUE)
                 )
    )
    pushViewport(
        viewport(name="viewport",
                layout.pos.row=1,
                layout.pos.col=1,
                xscale=range.x,
                yscale=range.y,
                clip=TRUE)
    )
    index <- 1

    gp <- gpar(col <- 3)

    # Draw lines for shape files.
    for(i in shp@polygons) {
        for(j in i@Polygons){
            name <- i@ID
            x <- j@coords[,1]
            y <- j@coords[,2]
            grid.polygon(x, y, default.unit="native",name=name, gp=gpar(fill="red"))
            index <- index + 1
        }
    }

    # Draw bounding box.
    name <- "bbox"
    x <- c(shp@bbox[1,1], shp@bbox[1,1], shp@bbox[1,2], shp@bbox[1,2], shp@bbox[1,1])
    y <- c(shp@bbox[2,1], shp@bbox[2,2], shp@bbox[2,2], shp@bbox[2,1], shp@bbox[2,1])
    grid.polyline(x, y, default.unit="native",name=name)
    upViewport(2)
}

自分のための「便利な自作関数」

このページのソースコードは、私がブログを書くために作ったもの。

使いたい人は勝手に使っても良いけれど、使用は自己責任で。バグがあるかも。

サポートに関しては、その時の気分次第。時間と心とお金に余裕があれば…。
ということで、返事が無くても怒らないように。疑問はコメント欄に。

今回のサンプルデータは以下の場所に。
https://docs.google.com/folder/d/0B9OtIs5BjRVhWm8zV1hCUFh6YnM/edit

以下からRの実行コード。

---

# Read a shape file.
X.geo <- readShapePoly('/home/yufujimoto/Desktop/japan.shp')

par(mfrow=c(1, 3),mar=c(2, 1, 4, 1))
cis.choropleth(X.geo, zcol="OCU", formula="scott", main="Number of Occupied Buildings")
cis.choropleth(X.geo, zcol="UNO", formula="scott", main="Number of Unoccupied Buildings")
cis.choropleth(X.geo, zcol="PRI", formula="scott", main="Average Buildings Rental Prices (Monthly)")

以下がRのソースコード。

---

# ======================================================
# cis.choropleth
# ------------------------------------------------------
# Description:
#   Make a choropleth map.
#
# Usage:
#   cis.colours(x, zcol, col1="orange", col2="red", formula="sturges", main="", cex=0.8)
#
# Arguments:
#   x       : SpatialPolygonsDataFrame
#   zcol    : A column used for making color ramp.
#   col1    : The name of color for lower bound.
#   col2    : The name of color for upper bound.
#   formula : Name of formula to make breakpoints.
#   main    : The title of this choropleth map.
#   cex     : Font size for the legend.
# ======================================================
cis.choropleth <- function(x, zcol, col1="orange", col2="red", formula="sturges", main="", cex=0.8, unit="", lty=1){
  # Read a required libraries.
  require(stats)
  require(grDevices)
  require(maptools)
  require(RColorBrewer)

  data <- x@data[which(colnames(x@data)==zcol)][,1]
  breaks <- cis.breaks(data, formula=formula)

  n <- breaks$class
  b <- breaks$breaks
  c <- cis.colours(col1=col1, col2=col2, class=n)
  l <- formatC(b, format="d")
  l.t <- "Legend"
  if(unit != ""){l.t <- paste(l.t,unit,sep=":")}

  plot(x,col=c[findInterval(data, b, all.inside=TRUE)],asp=1, lty=lty)
  grid()
  title(main=main)
  legend("topleft", legend=l, fill=c, cex=cex, bty="n", title=l.t)
}

# ======================================================
# cis.colours
# ------------------------------------------------------
# Description:
#   Return a list of interpolated colours with two colors
#   and with number of breakspoints.
#
# Usage:
#   cis.colours(col1="orange", col2="red", class=3)
#
# Arguments:
#   col1  : The name of color for lower bound.
#   col2  : The name of color for upper bound.
#   class : Number of breaks.
# 
# Value:
#   colours : A character with elements of 7 characters,
#             "#" followed by the red, blue, green. See 
#             'rgb' for more details.
# ======================================================
cis.colours <- function(col1="orange", col2="red", class=3){
  # Read a required library.
  require(grDevices)

  # r ramp from two colors.
  cols <- rgb(colorRamp(c(col1,col2))((0:class)/class),max=255)

  # Return a vector of colors.
  return(cols)
}

# ======================================================
# cis.breaks
# ------------------------------------------------------
# Description:
#   Get a vector of breaks with specifed formula.
#
# Usage:
#   cis.breaks(x, formula="sturges")
#
# Arguments:
#   x      : A data vector
#   formula: Name of formula to make breakpoints.
#
# Details:
#   .. $ formula: The name of used formula to make breaks.
#   .. $ class  : The suggested number of classes.
#   .. $ breaks : A vector of breakpoints.
# ======================================================
cis.breaks <- function(x, formula="sturges"){
  # Read a required library.
  require(stats)

  # Define a list to return.
  res <- list(formula=formula, class=0, breaks=0)

  # Get a name of formula.
  res$formula <- formula
  
  # Switch formula for making breaks.
  res$class <- switch(formula,
      sturges = nclass.Sturges(x),
      scott   = nclass.scott(x),
      FD      = nclass.FD(x)
  )

  # Get a vector of breaks with specifed formula.
  res$breaks <- quantile(x, probs=seq(0,1,1/(res$class)))

  # Return the result.
  return(res)
}