地図、地球儀

国名と形式を選んで地図データがダウンロードできます。形式には R があります。R形式は .RData（Rの作業スペース）です。Rで「ファイル」「作業スペースの読み込み」で読み込めます。
iraq を選んでダウンロードして読み込んで見ると単一のオブジェクト gadm が含まれます。plot(gadm)で表示できます。ダウンロードできるファイルはレベル0、１、２と3種あり、レベル0は国境線のみでした。後は行政区画の細かさのようです。
河川や湖沼は描かれません。
http://www.gadm.org/

.RData 以外の形式は、Shapefile、 ESRI file geodatabase、ESRI personal geodatabase、Google Earth KML があります。
ESRI 社の ArcGIS や ArcMAP が実質的な業界標準のようです。
QGISではShapefile、ESRI personal geodatabase を「レイヤの追加」で読み込むことが出来ます。
ESRI personal geodatabase は単一の .mdb ファイルです。
Shapefile、 ESRI file geodatabaseは複数のファイルから構成されています。

Gadmには国別以外に「Whole World」があります。Shapefile、 ESRI file geodatabase、ESRI personal geodatabase の解凍した後の大きさは、それぞれ 1.4GB、608MB、1.19GBでした。
R言語には、rgdal、PBSmapping、sp、maptools、gpclib と言ったパッケージがあり、これらのデータが扱えることになっています。
しかし、一度にデータフレームに読み込むもので、2GB程度のパソコンでは不可能でした。
ESRI personal geodatabase は .mdb なのでAccess 等で内容が見られます。データは海岸線や国境、行政区画などのように階層をなしていないようです。データは最小の行政区画のポリゴンの集合です。これを国名や国、地域コードで引くことで、国別や地域の地図を描きます。ポリゴンにはいくつかの属性が付与されていますが図形を描く際のレイヤの替わりになるものはないようです。ENGTYPE_1には Province や State と記されていますが、各国には1通りしかないので英語で何と呼ばれるかの参照にしかなりません。州名は有意義な情報で中国語表記もあります。これらの属性は元の情報にあったものを残したもので統一的な基準で含まれている訳ではないものと思います。

R言語には地図情報を扱ういくつかのパッケージがあります。
sp、maptools が基本的なパッケージのようです。shapefile形式でIRQ_adm.zipをダウンロードして試してみました。

shapefile を読み込む関数は複数あり、線画、ポリゴンなどで読み込むようです。

library(maptools)
#a <- readShapeLine("IRQ_adm/IRQ_adm2.shp")
#a <- readShapePoly("IRQ_adm2.shp")
a <- readShapeSpatial("IRQ_adm/IRQ_adm2.shp")
par(mar=c(0,0,0,0))
plot(a) # 地図を描く

国境線や海岸線が識別できる訳ではないので、国境を書くにはポリゴンをマージすることになります。gpclib にはその機能があるようです。ポリゴンデータを "gpc.poly" と言う型にして union() を呼び出します。

library(gpclib)
LM <- a@bbox # 図形の範囲
pos <- slot(a,"polygons") # データフレーム中の図形部分
n <- length(pos) # 図形データの総数
u <- as(pos[[1]]@Polygons[[1]]@coords,"gpc.poly")
for(i in 1:n)
{
po <- as(pos[[i]]@Polygons[[1]]@coords,"gpc.poly")
u <- union(u,po)
}
par(mar=c(0,0,0,0))
plot(u)

ポリゴンには11の属性が与えられています。しかし、有意そうなのは行政区画の名前だけです。行政区画の名前を含めて、属性はデータ作成時に元データから与えられたものだと思います。国が替われば属性の意味合いも変わるもので一律の基準で出来ているわけではないようです。
行政区画で選択して部分的な図が描けます。

par(mar=c(0,0,0,0))
plot(0,0,col="white",xlim=LM["x",],ylim=LM["y",],
xlab="",ylab="",axes=F)
for(i in which(a@data$NAME_2=="Basrah" |
a@data$NAME_2=="Abu al Khasib" |
a@data$NAME_2=="Shatt Al Arab" |
a@data$NAME_2=="Chibayish" |
a@data$NAME_2=="As Salman" |
a@data$NAME_2=="Al Zubair"))
{
polygon(pos[[i]]@Polygons[[1]]@coords)
}

Whole World データ

前述のように Global Administrative Areas は 2GB の実装のパソコンでは開けませんでいた。QGISなどでは表示するだけはなとか出来ますが操作は実用範囲外です。しかしメモリを足すだけで自由に地図が使えるならありがたいことです。
地図データとしては読めませんがデータベース（MDB）としては読めるので、地球を表すデータはどのようなものか gadm2.mdb を調べて置きます。（最大のテーブル gadm2 はフィールドを制限して部分的に読むことができますが全フィールド（*）を使うとR言語ごと落ちます。）

まず、データは8バイトの浮動小数点数（double）です。地球楕円体（準拠楕円体）の赤道半径は 6,378,137 m です。メートルのスケールで、赤道の一周は8桁の数字です。したがって、4バイト(float)では不足なことは確かです。しかし、メートルのオーダーで見た現実の地球は常に変化しているのだと思います。また、描画精度の面でもdoubleの制度は冗長です。
型を気にしないで済む点で double型は合理的な選択ですがプログラミング上は工夫の余地があります。R言語の演算精度は16桁程度のようです。

gadm2.mdb の図形データの行数は 218,238 です。行は行政区画を表すもので、国名や行政区画名が振られています。行は複数のポリゴンで出来ています。行の図形データのサイズは112から9,007,976byteです。
地球上の行政区画の総数は22万ほどのようです。
行政区画には以下のような区分が設定されています。各国で付与した区分の英訳だと思います。国が異なれば意味合いも変わるものだと思います。

> library(RODBC)
> mdb <- odbcConnectAccess("gadm2.mdb",readOnlyOptimize=T)
> a <- sqlQuery(mdb,"select DISTINCT ENGTYPE_1 from gadm2")
> a[[1]]
[1] <NA>
[3] Autonomous Republic Administrative Area
:

gadm2 テーブルのOBJECTIDとShapeフィールドを読んで図形データのサイズや行数を調べます。

> nrow <- sqlQuery(mdb,"select COUNT(*) from gadm2")
> nrow <- as.numeric(nrow)
> odbcQuery(mdb,"select OBJECTID,Shape from gadm2",rows_at_time =1)
[1] 1
> oid <- rep(NA,nrow)# 行のOBJECTID
> ln <- rep(NA,nrow) # 各行の図形領域のサイズ
> tp <- 0 # 緯度経度の総数
> ts <- 0 # ポリゴンの総数
> for(i in 1:nrow)
+ {
+ a <- odbcFetchRows(mdb, max = 1, buffsize = 32000)
+ oid[i] <- as.character(a[[1]][1])
+ u <- unlist(a[[1]][2])
+ ln[i] <- length(u)
+ tp <- tp + readBin(u[41:44],"int")
+ ts <- ts + readBin(u[37:40],"int")
+ }
> ts # ポリゴンの総数
[1] 315565
> tp # 緯度経度の総数
[1] 61316967
> range(ln) # 図形データの最小最大
[1] 112 9007976
> oid[which.max(ln)] # 最長の図形データの行のOBJECTID
[1] "34417"

地表の [経度、緯度] が示されている箇所は 6,132万地点あり、32万のポリゴンで表されています。

最も複雑な行政区画は OBJECTID=34417の行で、9MBの図形データで表されています。
ほとんどの行の図形データが1MB以下で、4MBを超えるものは他にないようです。

> odbcQuery(mdb,paste("select OBJECTID,NAME_0,NAME_1",
+ "from gadm2 where OBJECTID=34417"),
+ rows_at_time =1)
[1] 1
> a <- odbcFetchRows(mdb, max = 1, buffsize = 32000)
> str(a)
List of 2
$ data:List of 3
..$ : int 34417
..$ : chr "Chile"
..$ : chr "Magallanes y Antartica Chilena"

この場所はチリのマガジャネス・イ・デ・ラ・アンタルティカ・チレーナ州のようです。上部左はウェリントン島のようです。多数の島が詳細に描かれています。

DrawPoly <- function(shape)
{
u <- unlist(shape)
parts <- readBin(u[37:40],"int")
points <- readBin(u[41:44],"int")
# パートはdoube単位のオフセットだがポイント数として使用
pp <- rep(NA,parts)
p <- 45
for(i in 1:parts)
{
pp[i] <- readBin(u[p:(p+3)],"int")
p<- p + 4
}
if( parts > 1 )
pp[1:(parts-1)] <- pp[2:parts]-pp[1:(parts-1)]
pp[parts] <- points - pp[parts]
# パート（ポリゴン）を描画
for(k in 1:parts)
{
x <- rep(NA,pp[k])
y <- rep(NA,pp[k])
for(i in 1:pp[k])
{
x[i] <- readBin(u[p:(p+7)],"double")
y[i] <- readBin(u[(p+8):(p+15)],"double")
p <- p + 16
}
polygon(x,y)
}
}
#------------------------------------------------------
# 最も複雑な行政区画を描画
odbcQuery(mdb,"select Shape from gadm2 where OBJECTID=34417",
rows_at_time =1)
a <- odbcFetchRows(mdb, max = 1, buffsize = 32000)
par(mar=c(0,0,0,0))
plot(0,0,col="white",xlim=c(-75.74514,-72.50806),
ylim=c(-52.75875,-48.59431))
DrawPoly(a[[1]][1])

日本の図形データは1813行分あります。各行には日本語で行政区画の名称が付けられています。

mapdataパッケージ

R言語のパッケージは座標などから自動的データを取ってきて表示するような発想のようです。
mapdataパッケージがどのようにデータを持っているのか分かりませんが簡便に地図を描いてくれます。
河川については World Rivers Map Database と言うデータがあるようです。

> nrow <- sqlQuery(mdb,paste("select COUNT(*) from gadm2",
+ "where NAME_0='japan'"))
> nrow
Expr1000
1 1813
> d <- sqlQuery(mdb,paste("select NAME_0,NL_NAME_1,NL_NAME_2",
+ "from gadm2",
+ "where NAME_0='japan'"))
> str(d)
'data.frame': 1813 obs. of 3 variables:
$ NAME_0 : Factor w/ 1 level "Japan": 1 1 1 1 1 1 1 1 ...
$ NL_NAME_1: Factor w/ 47 levels "愛知県","愛媛県",..: 44 ...
$ NL_NAME_2: Factor w/ 1756 levels "あきる野市","あさぎり町",..

日本は [経度、緯度] が 870,940 記録され、3,225ポリゴンで描かれるようです。

nrow <- sqlQuery(mdb,paste("select COUNT(*) from gadm2",
"where NAME_0='japan'"))
nrow <- nrow[[1]]
odbcQuery(mdb,paste("select Shape from gadm2",
"where NAME_0='japan'"),
rows_at_time =1)
ln <- rep(NA,nrow) # 各行の図形領域のサイズ
tp <- 0 # 緯度経度の総数
ts <- 0 # ポリゴンの総数
par(mar=c(0,0,0,0))
plot(0,0,col="white",xlim=c(122.93320,153.98694),
ylim=c(24.04542,45.52279))
for(k in 1:nrow)
{
a <- odbcFetchRows(mdb, max = 1, buffsize = 32000)
DrawPoly(a[[1]][1])
u <- unlist(a[[1]][1])
ln[k] <- length(u)
tp <- tp + readBin(u[41:44],"int")
ts <- ts + readBin(u[37:40],"int")
}
ts # ポリゴンの総数
tp # 緯度経度の総数
range(ln) # 図形データの最小最大

日本の範囲は [北緯24.04542°、東経122.93320°] と [北緯45.52279°,東経153.98694°] によって表される矩形領域の中にあるようです。
愛知県をキーにして描くと左図下のようになります。

par(mar=c(0,0,0,0))
plot(0,0,col="white",xlim=c(136.708,137.7249),
ylim=c(34.6,35.42088))
for(id in objid)
{
odbcQuery(mdb,paste("select Shape from gadm2",
"where OBJECTID=",id),
rows_at_time =1)
a <- odbcFetchRows(mdb, max = 1, buffsize = 32000)
DrawPoly(a[[1]][1])
}

スペースシャトル地形データ

Shuttle Radar Topography Mission (SRTM) はNASAが行ったスペースシャトルからの測定に基いた標高データのようです。
N35E138.hgt.zip のように緯度と経度を名前に持った多数のファイルから出来ています。精度が1秒と3秒のデータがあるようです。無償で利用できるのは3秒のもののようです。SRTM1 や STRM3と区別されているようです。また、SRTM 90m と言う表記もあります。これは、緯度の3秒、赤道付近の経度の3秒が90m程度であることから来ているものと思います。

> library(raster)
> rs <- raster("N35E138.hgt/N35E138.hgt")
> a <- readAll(rs)
> plot(a)
> range(a@data@values,na.rm=T)
[1] -16 3742
> length(a@data@values)
[1] 1442401
> a@data@values[1:10]
[1] 963 947 911 869 813 782 757 751 736 742
> fi<-file.info("N35E138.hgt/N35E138.hgt")
> fi$size
[1] 2884802
> 1201*1201*2
[1] 2884802

N35E138.hgt.zip を解凍すると、N35E138.hgt になります。ファイルのサイズは2,884,802バイトです。これは、1201 x 1201 x 2 = 2,884,802 のようです。
1° を3秒で割ると1200です。データは2バイトの整数です。
R言語の raster パッケージは、このファイルを読み込んで図のような表示をします。
N35E138は、N35-36、E138-139 のデータのようです。
データの長さは 1,442,401 でファイル全体が標高データとして参照されます。ファイルの先頭を確認すると、Big Endian で読み込まれたデータと同じです。
データの最大は3742です。これは富士山の標高に近く3742 m と解釈して良さそうです。
また、データ中には 2417 の NA があります。どんな箇所が NA なのかは分かりません。ファイル上は16ビット整数なので NA は表現できませんが、0x8000 をNAとしているようです。16ビット整数の最小値の -32768 です。

library(raster)
rs <- raster("N35E138.hgt/N35E138.hgt")
a <- readAll(rs)
x <- rep(1:1201,1201) # x: 1:1201の数列を1201繰り返す
# y: 1を1201個、2を1201個、・・・、1021を1201個
y <- NULL; for(n in 1:1201) y <- c(y,rep(n,1201))
# 標高データの NA と負をゼロに
z <- ifelse(is.na(a@data@values), 0, a@data@values)
z <- ifelse(z < 0, 0, z)
# 標高データを 1:10 に
z <- floor(z / max(z)*9) + 1
z <- rainbow(10)[z] # 標高を色に
par(mar=c(3,3,2,1),mgp=c(1.5,0.5,0))
plot(x,y,col=z) # 描画

1列に並んだ数字列の順番の解釈を確認します。XY平面に、(0,0) から (0,1201)、(1,0) から (1,1201)、・・・、(1201,0) から (1201,1201)と描画して見ます。重い処理ですが1度きりなので良しとします。
結果は plot(a) でパッケージが描くものと異なっています。差異は単にY方向が逆なだけのようです。
したがって、N35E138.hgt ファイルのデータは、最も北から、「西から東へ」を繰り返して格納されているようです。表示のイメージそのままと理解出来ます。

北緯の最大は60°までのようです。N60E163.hgt を見るとデータの最初の方はNAが詰められています。実際に有効なのは60.3°台までのようです。データ数は緯度に関係なく一定のようです。
1°の表す東西の長さは緯度で大きく変わります。赤道では111km、60°では56km 程度です。緯度の絶対値が大きくなるほど水平方向のデータの密度は高いことになります。

標高や海抜には、１）最寄の港湾の平均海水面、２）各国の定めるジオイド、３）世界規模のジオイド、４）衛星軌道に基いた楕円体高、と言った異なる基準があることになります。
このデータはEGM96ジオイド面からの高さに加工したもののようです。
富士山の標高は 3,776 m とされています。データ中の最大値は 3,742 m でした。34m 足りませんがピークが測定されなくても不思議はありません。また、この位置のジオイド高は 40m ほどで、小さすぎるので楕円体高ではないと考えられます。
富士山の標高は日本で定めたジオイド面から測られたものですが、EGM96との差は 1.18 m 程度のようです。EGM96の方が僅かに標高が高く計算されることになります。

標高について

国内で使用されている山の高さや標高は東京湾の平均海水面からの高さです。これは世界規模の規準ではないと考えられます。

日本の標高もジオイド面からの高さである点では世界規模の規準と変わりありません。ジオイドのゼロを東京湾の平均海水面としている点で異なっています。
現在使用されてている高さは、楕円体高と標高の2種類のようです。衛星を使用して測定される楕円体高は数mm 程度の精度があるようです。一方、標高はジオイド面からの高さで自明ではありません。
日本ではジオイド面の規定を「日本のジオイド」や「日本のジオイド2000」と言うようです。これは「JGEOID2000」との差として与えられているようです。「JGEOID2000」は世界規模のジオイド高データである「EGM96」を日本部分について詳細化にしたもののようです。「JGEOID2000」は「ジオイド高」のデータで、「日本のジオイド2000」は「東京湾の平均海水面規準への修正量」のデータのようです。この修正量の最小値は -0.384 m、最大値は0. 358 m のようです。、「日本のジオイド2000」は、楕円体高の測定値から水準点の標高を引いて算出されたもののようです。「GSIGEO2011+2000」と言うジオイド高のデータが2kmメッシュで提供されているようです。

国土地理院のつくば市の「地球ひろば」には展示用の三角点（標高28.0m）があり、正確な緯度、経度が示されています。

EGM96 の計算値より 1m18cm ほど大きく、前述の40cm弱の差ではないようです。この地点の楕円体高は 66.91m から 68.10m ほどと言うことになります。
前述のSRTMデータの N36E140.hgt の、この場所の高さは 25 m です。位置は完全に一致するものではないので、周囲の平均を取ると 27.33 m です。EGM96の方が低い位置から測るので大きな値になると考えましたが実際には違うようです。

気象庁の富士山頂の地震計の設置場所は標高が 3772 m と記されています。

前述のようにSRTMの北緯35°、東経123°のデータの最大値は 3742 でした。
このことは、SRTM の標高は、楕円体高ではないことは間違いなさそうです。文書としては単位がメートルであり、「 WGS84/EGM96 geoid」を参照とあります。

World Rivers Map Database

河川は地図に不可欠に思えます。R言語の mapdata パッケージには World Rivers Map Database と言う説明がありました。
このパッケージの map() は、最初の引数がデータベース名になっていて、そのデータベースの全体、または部分を表示するもののようです。

> library(mapdata)
> d <- map("world")
> str(d)
List of 4
$ x : num [1:15356] -130 -132 -132 -133 NA ...
$ y : num [1:15356] 55.9 56.7 57 58.4 NA ...
$ range: num [1:4] -180 190.3 -85.4 83.6
$ names: chr [1:2284] "Canada" "South Africa" "Denmark" "Great Lakes:Superior, Huron, Michigan" ...
- attr(*, "class")= chr "map"
> d <- map("rivers")
> str(d)
List of 4
$ x : num [1:1357] -123 -122 -123 -125 -118 ...
$ y : num [1:1357] 65.2 66.2 66.2 66.2 66.4 ...
$ range: num [1:4] -164.5 160.8 -36.9 72.9
$ names: chr [1:239] "1" "2" "3" "4" ...
- attr(*, "class")= chr "map"

"world" と言うデータベースには2284の国や地域があり、"Canada"、"South Africa" と言った国名で表示範囲を選択できます。
"rivers" と言うデータベースには 239 の河川の図形があり、番号で識別されているようです。
全体を描くと以下のようになります。

map("world",mar=c(0,0,0,0),mgp=c(0,0,0),
col="blue",interior=F)
map("rivers",add=T,col="red")

チグリス、ユーフラテス川流域を表示する場合は以下のようにします。

map("world",col="blue",
region=c("Iraq","Iran","Syria",
"Kuwait","Saudi Arabia"))
map("rivers",add=T,col="red")

世界地図の全図と、それと重なる河川の大きな図があって、その一部分だけが見えているイメージのようです。
国境などの線は interior=F で非表示になります。これは、データが独立しているのではなく、前述の union() に依るもののようです。
河川データは値から [経度、緯度] のようです。

データベースには他に、japan、china、nz、world2 があります。
日本の全図は、worldのjapan とjapan の2通りあることになります。

map("world","japan")
map("japan")

map("japan","Aichi")

riverse には日本の河川データはありません。
国土地理院は河川データを公開しているので利用する方法はあるのだと思います。

この世界全図は 15,356 ポイントとコンパクトです。このデータで地球儀を描いて見ます。 rgl パッケージで地図を描けば平面へのプロジェクションの問題を回避できます。個人がパソコンで見る分には平面にする必要性は特にないと思います。
rgl でプロットした図は自由に回転したり拡大縮小出来ます。

このデータには 3902 ポイントの NA があります。世界地図は一筆書きできるものではないので複数の線データをNAで区切っているものと思います。したがって、クローズしている線は始点と終点が同じ値になっていることになります。この閉じた線はポリゴンで塗ることが出来ると思います。
また、経度は東経を正、西経を負で表していますが東経には180°を超えるものがあります。189.8° ほどまであり、連続した線を描くために調整されているのだと思います。

RGLのXYZ直交座標系で描くために [経度、緯度] を（x、y、z）に変換します。緯度は「地理緯度」だと解釈しました。準拠楕円体上の接線の法線と、赤道面のなす角です。地心座標に変換して、地心を通る直線と楕円体表面の交点の座標を計算します。

# 度で与えられた経度と緯度（地理緯度）からXYZ座標を計算
XYZ <- function(la,fa) #!!! 経度、緯度の順
{
a <- 6378137 # WSG84準拠楕円体の赤道半径
b <- 6356752.3142 # f:298.257223563から極半径
a2 <- a^2
b2 <- b^2
ph <- atan(b2/a2*tan(fa*pi/180)) # 地理緯度から地心緯度に
la <- la * pi / 180 # 経度を弧度に
# tan(ph)の傾きの直線と楕円の交点
ta <- tan(ph)
ta2 <- ta^2
u <- sqrt(b2 + a2*ta2)
xx <- a*b / u
c(xx*cos(la),xx*sin(la),a*b*ta / u)
}

library(rgl)
library(mapdata)
d <- map("world",plot=F)
n <- length(d$x)
xyz <- array(NA,c(n,3))
for(i in 1:n)
{
if( !is.na(d$x[i]) )
{
xyz[i,] <- XYZ(d$x[i], d$y[i])
}
}
plot3d(xyz[,1],xyz[,2],xyz[,3],t="l",col="blue",axes=F)
# 赤道
a <- 6378137
t <- seq(0,2*pi,length.out=100)
plot3d(a*cos(t),a*sin(t),0,t="l",col="red",add=T)

マウス操作で回転して北海道の辺りを拡大すると図のように見えます。

rgl には polygon3d() があって塗れるかもしれないと考えましたが簡単ではないようです。図形データには線が交差しているものが含まれています。
準拠楕円体を塗ることを考えます。

楕円体の図形を生成する関数は cdaパッケージにありました。
上のスクリプトの後半の描画部分を以下に差し替えると左図になります。

library(cda)
a <- 6378137
b <- 6356752.3142
open3d()
shapelist3d(rgl.ellipsoid(a=a,b=a,c=b),
col="yellow")
plot3d(xyz[,1],xyz[,2],xyz[,3],t="l",
col="blue",axes=F,add=T)
# 赤道
t <- seq(0,2*pi,length.out=100)
plot3d(a*cos(t),a*sin(t),0,t="l",
col="red",add=T)

山脈

山脈も河川と共に歴史の話に良く登場します。ヒマラヤは、ほぼ完全に人の往来を妨げていたのだと思います。また、コーカサスやアルプス、ピレネーは軍事行動を制限していました。前8世紀にはウクライナのキンメリア人を追ったスキタイがコーカサスの東端を通ってメソポタミア（メディア）へ侵入したと考えられているようです。BC216、ハンニバルは象を伴ってアルプスを越えローマに迫りました。

人口の大半が標高の低いところにあることは確かそうです。平らな広い場所は、低地なら平野、標高が高いところは高原と呼ばれるようです。デカン高原の標高は 300-600m のようです。
おそらく、数百メートルを越える標高の場所は人の移動に障害となるのだと思います。
最も標高の高い都市はボリビアの都市で、エル・アルト、ポトシが上げられています。エル・アルトは標高4,150mで人口は83万に程度のようです。

SRTM30 のデータは、南緯60° から北緯90° で、経度は360°分揃っているようです。1つのファイルは、緯度50°、経度40°分です。3 x 9 = 27 ファイルから成ります。ファイル名は、「西端北端.dem」となっています。e020n40.dem は、東経20°から東経60°、南緯10°から北緯40°の区域の標高です。w020n40.dem は、西経20°（-20°）から東経20°の範囲になります。
ドキュメントでは、このデータの最小値は -415 m、最大値は 8685 m です。
データの格納位置から計算した最小値の [経度、緯度 ] は、[35.50833、31.31156] です。死海（東経35度30分0秒、北緯31度30分0秒）のようです。
最大値の [経度、緯度 ] は、[86.84167、27.975] です。エベレスト（東経86度56分40秒、北緯27度59分16秒）の近傍と言うことです。

人の移動を妨げてきた山脈を地図に加えたいと思います。地球規模で考えて、2,000 m 以上を描いてみようと思います。

R言語の contour() で、SRTM30の e060n40 を描くと左図のようになります。
期待したのは2000m以上の地帯を囲むポリゴンですが細かな図形の集まりになってしまいます。
データ値が2000m以上を赤、0m以下を青に塗ると下図のようになります。欲しいのは赤の部分の輪郭のポリゴンです。それが得られれば、polygon3d()を使って、地球儀の上に山脈が描けそうです。
しかし試行錯誤して見て単純な話しではないことを知りました。

まず、等高線に相当する図形データを得ることが必要です。contour() は図形データを返さないようなので、等高線に相当する図形を描くR言語のスクリプトを作ってみました。一般に等高線は、指定の高さの (x,y) を計算します。これに対して今回は格子状の (x,y) のどこを通るかと言う考え方で作ります。 (x,y) は [経度、緯度] で、精度面でも差分を取って(x,y) を調整する必要性はないと考えました。

# 探索方向による確認順
D <- matrix(c(
-1,-1, -1, 0, -1, 1, 0, 1, 1, 1, 1, 0, 1,-1, 0,-1,
-1, 0, -1, 1, 0, 1, 1, 1, 1, 0, 1,-1, 0,-1, -1,-1,
-1, 1, 0, 1, 1, 1, 1, 0, 1,-1, 0,-1, -1,-1, -1, 0,
0, 1, 1, 1, 1, 0, 1,-1, 0,-1, -1,-1, -1, 0, -1, 1,
1, 1, 1, 0, 1,-1, 0,-1, -1,-1, -1, 0, -1, 1, 0, 1,
1, 0, 1,-1, 0,-1, -1,-1, -1, 0, -1, 1, 0, 1, 1, 1,
1,-1, 0,-1, -1,-1, -1, 0, -1, 1, 0, 1, 1, 1, 1, 0,
0,-1, -1,-1, -1, 0, -1, 1, 0, 1, 1, 1, 1, 0, 1,-1),
ncol=16, nrow=8, byrow=T)
# 検出ポイントによる探索方向
E <- matrix(c(6,5,4,
7,0,3,
8,1,2), ncol=3, byrow=T)
# 接続ポイントの探索
FindNext <- function(cs, rs, c, r, dr)
{
for(i in seq(1,by=2,length.out=8))
{
tc <- c + D[dr,i]; tr <- r + D[dr,i+1]
nr <- which(rs==tr)
if(length(nr)>0)
{
nc <- which(cs[nr]==tc)
if(length(nc)>0)
{
return(c(E[D[dr,i]+2,2-D[dr,i+1]],nr[nc]))
}
}
}
return(NULL)
}
#-------------------------------------------------
# 等高線（境界のポイントを結ぶ）
SimpleContour <- function(m,level)
{
ncol <- dim(m)[1]
nrow <- dim(m)[2]
# 変化点の抽出
m[which(m<level)] <- 0
m[which(m!=0)] <- 1
# 行、列の両方で境界を検出するほうが良いが容量不足になる
m[,1:(nrow-1)] <- xor(m[,1:(nrow-1)],m[,2:nrow])
#m2 <- m
#m2[,1:(nrow-1)] <- xor(m[,1:(nrow-1)],m[,2:nrow])
#m[1:(ncol-1),] <- m2[1:(ncol-1),]
# + xor(m[1:(ncol-1),],m[2:ncol,])
ps <- which(m>=1,arr.ind=T)
rm(m)
cs <- ps[,1]
rs <- ps[,2]
#
plot(0,0,col="white",xlim=c(0,ncol),ylim=c(0,nrow))
#plot(cs,rs,col="green",xlim=c(0,ncol),ylim=c(0,nrow))
#
x <- rep(NA,ncol*6)
y <- rep(NA,ncol*6)
dr <- 1 # 1下,2,3,4,5,6,7,8(反時計周り)
while(length(cs)>1)
{
c <- cs[1]; r <- rs[1]
cs[1] <- 0; rs[1] <- 0
x[1] <- c; y[1] <- r
p <- 1
dr <- 1
repeat
{
n <- FindNext(cs, rs, c, r, dr)
if( length(n) == 0 ) break
dr <- n[1]
c <- cs[n[2]]; r <- rs[n[2]]
cs[n[2]] <- 0; rs[n[2]] <- 0
p <- p + 1
x[p] <- c; y[p] <- r
}
if(p>2)
{
k1 <- abs(x[p]-x[1]); k2 <- abs(y[p]-y[1])
if( (k1<2) && (k2<2) )
{
polygon(x[1:p],y[1:p])
} else {
lines(x[1:p],y[1:p])
}
}
n <- which(cs>0)
cs <- cs[n]
rs <- rs[n]
}
}
#-------------------------------------------------
ncol <-4800
nrow <- 6000
fn <- "srtm30/E060N40.DEM"
f <- file(fn, "rb")
m <- array(NA,c(ncol,nrow))
for(i in 1:nrow)
{
a <- readBin(f,integer(),ncol,size=2,endian="big")
m[,nrow-i+1] <- a
}
close(f)
rm(a)
windows(3,3)
par(mar=c(3,3,2,1),mgp=c(1.5,0.5,0))
SimpleContour(m,2000)
title("素朴なスクリプトで描画(SRTM30区画)",cex.main=0.8)

境界自体は排他的論理和や、水準を引いた後の2点の値の積の符号で検出できます。容量不足なので境界の座標データだけを残して元データは解放することにしました。輪郭抽出アルゴリズムのように逐次周辺の8個の高さを見ないで済む範囲で接続しました。

21 x 21 極小さな部分を見ると、contour() は綺麗に線を引きます。このスクリプトでは細切れになってしまいます。

しかし、全体を描けば大差はないようです。
細切れになったものを上手く繋ごうと試行錯誤していてやっと気が付いたのは図形を連結して行くと輪郭抽出アルゴリズムと同じようにホールのできる可能性が生じることです。ホールのある図形を1つの図形として扱う方法が無いので連結しても意味がありません。
細かいポリゴンをたくさん描くことは必須のようです。
残る問題は、標高の高い山脈の中心部分のポリゴンが生成されていないことです。

そこで、東西方向に標高の高いところの連続区間を取るようなアルゴリズムに変えました。

しかし、残念ながら polygon3d() で塗る目論見は上手く行きませんでした。
下図は polygon3d() ですが、fill＝FALSE として描画したものです。塗ろうとするとエラーが発生します。

まず、NAで区切って複数図形を一度に渡すことが出来ないようです。また、先頭と最後のポイントを一致させる必要はなく、逆にエラーになります。
これらを修正しても十数個に1つの割合でエラーになるポリゴンがあります。
エラーになる最初のポリゴンは線分上に点が乗ってしまうケースでした。
各点について、その点を端に持たない全ての線分に付いて、領域が重なり、線分の一端と結んだ線の傾きが同じものを探します。あれば、凹凸をなまらせる方向に1つ点を移動して調整しました。この操作は1つの直線上に3つ以上点があると複雑になるので図形の単純化と繰り返しで変化が無くなるまで繰り返しました。
それでも数十個に1つ程度でエラーになります。ポリゴンの形状に異常はありません。良く分かりませんが細長い狭い箇所のある図形がエラーになるようです。エラーは三角形に分割することに失敗すると言うことのようです。 polygon3d() は x、y、z を対等に扱っているのではなく、x、y だけで三角形に分割しているようです。したがって、XY平面に垂直な図形はエラーになりますが、ここでは該当しません。

もう一つの問題はポリゴンを単純化したこともあって、アルゴリズムから東西に長い線分を生じることです。直線は地表にもぐってしまうのでヒマラヤ山脈の中心部分が見えなくなってしまいました。改めて線分の長さが 1° を越えたら分割して描きました。
ここまで試して、3Dについてもう少し知る必要がありそうだと考えました。「rgl パッケージ」

3D Quadrangle Mesh

rgl パッケージの使い方を「rgl パッケージ」のように試して見ました。
qmesh3d(バーテックス、インデックス) として四角形メッシュ図形が作れます。
SRTM30のデータは、東西に40°、南北に50°を1ファイルとして作られています。このファイルの単位でバーテックスとインデクスを抽出し qmesh3d を作って連結しました。

2000mを越える地帯は少なく、ヒマラヤ山脈は例外のようです。中国との交通は大きく迂回する以外には不可能だったと考えて良さそうです。1万1千年前ころに間氷期に入りますが、ヒマラヤの周辺は長い間、氷河による水の豊かな場所だったのだろうと推測できます。
ウクライナにいたスキタイがメディアに侵入したのはコーカサス山脈を迂回するのにカスピ海沿岸を通ったことが理由のようです。
セレウコスはチャンドラグプタに東部地域を割譲しました。アショーカ王のギリシア文字・アラム文字碑文はカンダハールにあると言うことです。カンダハールは標高千メートル強のようです。
古代の都市を [経度、緯度] でプロットすれば河川や山脈との位置関係がはっきりします。

# 度で与えられた経度と緯度（地理緯度）からXYZ座標を計算
XYZ <- function(la,fa) #!!! 経度、緯度の順
{
if(is.na(la)||is.na(fa))
return(c(NA,NA,NA))
a <- 6378137 # WSG84準拠楕円体の赤道半径
b <- 6356752.3142 # f:298.257223563から極半径
a2 <- a^2
b2 <- b^2
ph <- atan(b2/a2*tan(fa*pi/180)) # 地理緯度から地心緯度に
la <- la * pi / 180 # 経度を弧度に
# tan(ph)の傾きの直線と楕円の交点
ta <- tan(ph)
ta2 <- ta^2
u <- sqrt(b2 + a2*ta2)
xx <- a*b / u
c(xx*cos(la),xx*sin(la),a*b*ta / u)
}
#------------------------------------------------------
# 色を塗った地球楕円体を描く
library(rgl)
library(cda)
a <- 6378137
b <- 6356752.3142
plot3d(rgl.ellipsoid(a=a,b=a,c=b),xlab="",ylab="",zlab="",
axes=F, col="khaki")
clear3d(type="lights"); light3d(specular="black")
# mapdataパッケージのworld,riversを3Dで表示
library(mapdata)
d <- map("world",plot=F)# 行政区画
d <- mapply(XYZ,d$x,d$y)# XYZ座標系に変換
plot3d(d[1,],d[2,],d[3,],
t="l",col="gray40",add=T)
d <- map("rivers",plot=F)# 河川
d <- mapply(XYZ,d$x,d$y) # XYZ座標系に変換
plot3d(d[1,],d[2,],d[3,],
t="l",col="blue",add=T)
#------------------------------------------------------
# 赤道
t <- seq(0,2*pi,length.out=100)
plot3d(a*cos(t),a*sin(t),0,t="l",col="red",add=T)
# 緯線 -60,-10,40,
t <- seq(0,360,length.out=100)
for(fa in c(-60,-10,40))
{
d <- mapply(XYZ,t,fa)
lines3d(d[1,],d[2,],d[3,],col="gray60")
}
# 経線
las <- c(-140,-100,-60,-20,20,60,100,140,180)
t <- seq(-90,90,length.out=50)
for(la in las)
{
d <- mapply(XYZ,la,t)
lines3d(d[1,],d[2,],d[3,],col="gray60")
}
#------------------------------------------------------
# 標高 2000m 以上を描画
load("srtm30/world2k.R.data")
shade3d(world2k,color="#808000")
#------------------------------------------------------
# 都市名を入れる例
Num <- function(c)
{
c( c[1]+c[2]/60+c[3]/3600, c[4]+c[5]/60+c[6]/3600 )
}
cities <- list(
list(Num(c(33,5,37, 44,34,50)),"Ctesiphon","クテシフォン",c(0,0.6)),
list(Num(c(33,5,40, 44,31,20)),"Seleucia","セレウキア",c(1.1,0.5)),
list(Num(c(36,21,34, 43,09,10)),"Nineveh","ニネヴェ",c(1,0.5)),
list(Num(c(32,32,11, 44,25,15)),"Babylon","バビロン",c(1,1)),
list(Num(c(33,20,00, 44,26,00)),"Baghdad","バグダート",c(0.1,-0.1)),
list(Num(c(30,57,45, 46,06,11)),"Ur","ウル",c(1,1.5)),
list(Num(c(31,19,20, 45,38,10)),"Uruk","ウルク",c(1,1)),
list(Num(c(35,27,24, 43,15,45)),"Assur","アッシュール",c(1,0.5)),
list(Num(c(32,11,26, 48,15,28)),"Susa","スサ",c(0.5,0.5)),
list(Num(c(29,56,04, 52,53,29)),"Persepolis","ペルセポリス",c(0,0.5)),
list(Num(c(30,00,42, 52,24,28)),"Anshan","アンシャン",c(1,0.5)),
list(c(33.7560,72.8291), "Gandhara","ガンダーラ",c(0,0.5)),
list(Num(c(36,46,23, 66,52,25)),"Bactra","バクトラ",c(0,0.5)),
list(Num(c(31,37,01, 65,43,01)),"Kandahar","カンダハール",c(0,0.5)),
list(c(25.611, 85.144),"Pataliputra","パータリプトラ",c(0,0.5)),
list(c(31.198, 29.9192),"Alexandria","アレクサンドリア",c(1,0.5)),
list(c(40.76, 22.519167),"Pella","ペラ",c(1,0.5)) )
for(i in 1:length(cities))
{
xyz <- XYZ(cities[[i]][[1]][2],cities[[i]][[1]][1])
points3d(xyz[1],xyz[2],xyz[3],color="red")
text3d(xyz[1],xyz[2],xyz[3]*1.03,cities[[i]][[2]],
color="black",adj=cities[[i]][[4]],cex=0.8,
font=2)
}

デフォルトでは大きなスペキュラが描かれ中心部分の表示が飛んでしまいます。ライトの設定は一端クリアすると上手く行くので、何かまだ知らない設定があるようです。
テキストも表示位置によって地面にもぐりこむので工夫が必要です。

2000m以上の箇所の着色は、world2k.R.data をロードして行っています。load() は退避した変数を復元するもので、world2k を復元します。world2k.R.data に save() した変数は、以下の手順で作りました。

C#でプログラムを書いて、SRTM30 のファイルから、ファイルごとに座標とインデクスのファイルを作りました。

class ColRow
{
public int col; // ラスターデータのカラム、あるいは経度
public int row; // ラスターデータの行、あるいは緯度
public ColRow(int col, int row)
{
this.col = col; this.row = row;
}
}
const int MESHMAX = 120; // 1度
const int MESHMIN = 12; // 6分
const int MESHLIM = MESHMIN * MESHMIN / 2;
int[,] z;// MESHMIN x MESHMIN のlevel以上のカウント
// 経度は0から4799で、40度の範囲
// 緯度は0から5999で、50度の範囲
List<ColRow> lafa = new List<ColRow>();
List<int> indices = new List<int>();
string fileName;
public SRTM30QMesh(string fileName, int ncol, int nrow, int level)
{
this.fileName = fileName;
BinaryReader br = new BinaryReader(
File.OpenRead(fileName));
z = new int[ncol / MESHMIN, nrow / MESHMIN];
for (int r = 0; r < z.GetLength(1); r++)
{
for (int i = 0; i < MESHMIN; i++)
{
// 1行読み込み
byte[] by = br.ReadBytes(ncol * 2);
// 1行のカウント
int p = 0;
for (int c = 0; c < z.GetLength(0); c++)
{
// 1行のMESHMIN個のカウント
for (int k = 0; k < MESHMIN; k++)
{
int v = (sbyte)by[p++];
v = v * 256 + by[p++];
if (v >= level) z[c, r]++;
}
}
}
}
br.Close();
// z[,]の処理
ColRow cr;
while ((cr = Find()) != null)
{
CreateMesh(cr);
}
// 結果をファイルへ出力
Write();
}
// 座標とインデクスをそれぞれファイルに出力
public void Write()
{
string fn = System.IO.Path.GetFileNameWithoutExtension(fileName);
TextWriter tw = File.CreateText(@"C:\CodeSamples\R_WORK\"+fn+"_coord.dat");
foreach (var d in lafa)
{
tw.WriteLine(d.col + " " + d.row);
}
tw.Close();
tw = File.CreateText(@"C:\CodeSamples\R_WORK\"+fn+"_index.dat");
for (int i = 0; i < indices.Count; i += 4)
{
tw.WriteLine(indices[i] + " " + indices[i + 1] + " " + indices[i + 2] + " " + indices[i + 3]);
}
tw.Close();
}
// z[,]全体から出力対象を探し、見つかったら処理済にセットする
ColRow Find()
{
for (int i = 0; i < z.GetLongLength(1); i++)
for (int k = 0; k < z.GetLongLength(0); k++)
if (z[k, i] > MESHLIM)
{
z[k, i] = -1;
return new ColRow(k, i);
}
return null;
}
// 矩形を検出して登録
void CreateMesh(ColRow cr)
{
int w = 1;
int c = cr.col + 1;
while ((w < MESHMIN) && (c < z.GetLength(0)))
{
if (z[c, cr.row] <= MESHLIM) break;
w++;
z[c, cr.row] = -1;
c++;
}
int h = 1;
int r = cr.row + 1;
while ((h < MESHMIN) && (r < z.GetLength(1)))
{
bool high = true;
for (int k = 0; k < w; k++)
{
if (z[cr.col + k, r] <= MESHLIM)
{
high = false;
break;
}
}
if (!high) break;
h++;
for (int k = 0; k < w; k++)
z[cr.col + k, r] = -1;
r++;
}
int la1 = cr.col * MESHMIN;
int la2 = la1 + w * MESHMIN;
int fa1 = cr.row * MESHMIN;
int fa2 = fa1 + h * MESHMIN;
SetQMesh(la1, fa1, la2, fa2);
}
//---------------------------------------------------
// lafaリストの要素間の比較用関数
ColRow cur_lafa;
bool compareLaFa(ColRow lafa)
{
if (lafa.col != cur_lafa.col) return false;
if (lafa.row != cur_lafa.row) return false;
return true;
}
//---------------------------------------------------
// 四角形を登録
void SetQMesh(int x1, int y1, int x2, int y2)
{
int i;
Action<int, int> Add = (Action<int, int>)((x, y) =>
{
cur_lafa = new ColRow(x, y);
i = lafa.FindIndex(compareLaFa);
if (i < 0)
{
i = lafa.Count; // 新出のポイント
lafa.Add(cur_lafa);
}
indices.Add(i);
});
Add(x1, y1);
Add(x2, y1);
Add(x2, y2);
Add(x1, y2);
}
//---------------------------------------------------
}

# メッシュデータの連結してworld2k.R.dataを作る
MeshCat <- function(mesh1,mesh2)
{
mesh <- mesh1
mesh$vb <- cbind(mesh1$vb,mesh2$vb)
mesh$ib <- cbind(mesh1$ib,mesh2$ib + dim(mesh1$vb)[2])
mesh
}
# SRTM30の区画ごとのデータをメッシュデータに
SRTM30Mesh <- function(fn)
{
bLa <- as.integer(substr(fn,2,4))
bFa <- as.integer(substr(fn,6,7))
if(toupper(substr(fn,1,1))=="W") bLa <- -bLa
if(toupper(substr(fn,5,5))=="S") bFa <- -bFa
lafa <- scan(paste(fn,"_coord.dat",sep=""))
indx <- scan(paste(fn,"_index.dat",sep=""))
if(length(lafa)==0)
return(NULL)
m <-matrix(1,ncol=length(lafa)/2,nrow=4)
p <- 1
for(i in 1:(length(lafa)/2))
{
xyz <- XYZ(bLa+lafa[p]/120.0,bFa-lafa[p+1]/120.0)
m[1,i] <- xyz[1]; m[2,i] <- xyz[2]; m[3,i] <- xyz[3]
p <- p + 2
}
i <- matrix(indx+1,nrow=4)
qmesh3d(m,i)
}
# メッシュデータファイルの名前
fns <- c(
"E020S10","E060S10","E100S10","E140S10","W020S10","W060S10","W100S10","W140S10","W180S10",
"E020N40","E060N40","E100N40","E140N40","W020N40","W060N40","W100N40","W140N40","W180N40",
"E020N90","E060N90","E100N90","E140N90","W020N90","W060N90","W100N90","W140N90","W180N90")
# 連結操作
mesh1 <- SRTM30Mesh(fns[1])
mesh2 <- SRTM30Mesh(fns[2])
world2k <- MeshCat(mesh1,mesh2)
for(fn in fns[3:length(fns)])
{
mesh <- SRTM30Mesh(fn)
if(is.null(mesh)) {
cat(fn,"\n")
} else {
world2k <- MeshCat(world2k,mesh)
}
}
# 描画
shade3d(world2k,color="#808000")
# 保存
save(world2k,file="world2k.R.data")
str(world2k)

save(world2k,file="world2k.R.data") で保存しました。このファイルは4バイト整数に件数を掛けて計算した容量より僅かに小さいようです。
ここで作成したデータは荒いもので、1° か 6分の四角形メッシュで描かれています。

mikeo_410