スポンサーサイト

上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。

簡単な色分け地図

$$

data <- read.table("original.dat",header=F)
sd0 <- sd(data[,2])
m0 <- mean(data[,2])

result <- ((data[,2]-m0)/4/sd0 + 0.3)*1.1
result[result[]<0] <- 0
result[result[]>1] <- 1
mycol <- rainbow(100)
prefcol <- mycol[result*100]
result2 <- t(col2rgb(prefcol))
write.table(result2,"result.dat",row.names=F,col.names=F,quote=F,sep="/")

cpt <- mycol[c(1,25,50,75,100)]
result2 <- t(col2rgb(cpt))
write.table(result2,"cpt.cpt",row.names=F,col.names=F,quote=F,sep="/")



#!/bin/bash
chars=( 0/255/240 255/15/0 0/255/117 0/255/148 255/77/0 255/46/0 66/255/0 255/229/0\
189/255/0 112/255/0 51/255/0 0/255/224 255/0/15 0/255/179 0/25/255 255/46/0\
255/15/0 204/255/0 0/255/56 0/255/133 255/153/0 5/255/0 255/138/0 255/46/0\
66/255/0 0/10/255 97/255/0 219/255/0 0/255/10 255/184/0 255/31/0 0/71/255\
0/163/255 204/255/0 255/0/122 82/255/0 204/255/0 255/245/0 255/0/122 0/255/25\
0/133/255 66/255/0 255/214/0 255/229/0 0/255/41 255/184/0 )
export k=1
psxy pref/pref$k.dat -Jm0.25 -R127/150/29/46 -M -V -W1 -Bg5a5 -G200/200/0 -P -K > data-test.ps

for test in "${chars[@]}"; do
psxy pref/pref$k.dat -Jm -R127/154/29/46 -M -V -W1 -G$test -P -O -K >> data-test.ps

echo $k
k=$(($k + 1))
done

psxy pref/pref$k.dat -Jm -R127/154/29/46 -M -V -W1 -G200/200/0 -P -O -K >> data-test.ps
psscale -D4.8/1.7/3/0.2 -Cdentists_stat/cpt.cpt -O >> data-test.ps

convert -density 300 data-test.ps data-test.png

exit


data-test2-mini.png


上記の地図は、公開データ「国土交通省 国土数値情報ダウンロードサービス 行政区域(面)H25年」を利用し作成した。
スポンサーサイト

コメントの投稿


非公開コメント

コメント

Profile

zoa

Author:zoa
Metropolis (1927)

Calendar
05 | 2017/06 | 07
- - - - 1 2 3
4 5 6 7 8 9 10
11 12 13 14 15 16 17
18 19 20 21 22 23 24
25 26 27 28 29 30 -
Labels
Admin
Previous Posts
Recent Comments
Recent Trackback
Archive
このページのトップへ
上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。