salvageon
結城浩さんのcodeIQの問題 https://codeiq.jp/ace/yuki_hiroshi/q1215
フィードバックをいただきました。両方正解でした。
Rを使っていた人は2名だけのようでしたが,ビッグデータの取り扱いは統計ツールを日常的に使う人にも関係は深いと思うので,もっとR使いがいてもいいのではないでしょうか。私はgmpというパッケージで超長整数を扱って計算しました。
DB1は単調増加とわかった段階で,微分方程式の数値解析のオイラー法みたいな計算をしましたが,解説ではバイナリサーチをしていました。この場合はオイラー法のほうがAPIをたたく回数は減ると思いますが…
解答集
http://togetter.com/li/758302
以下にanswer.txtをさらします。(画面からはみ出る部分に改行を追加)
V406435859539156181269150751031 V1101943557675920722238136981003 ENV: 最終的に「R」。途中経過で電卓,Excel, 少しだけpython POINT: DB1はキーが単調増加。DB2はキーを100bitの2進数で表し いちばん下の1の桁を上下反転させたものに残りの桁を加えたもの の順に並んでいる。 (キーQを割り切れる最大の2のべき乗を2^k(k>=0)とした時, 2^(99-k)+(Q-2^k)/2^(k+1)がindexになる) Rのソースと実行結果
- -
- -
- -
- -
- -
- -
- -
- -
実はRで解く前にwindowsの電卓(超長整数を扱えるようだ)を使っていろいろ試していたのですが,DB2では2進数と10進数の変換に以下のサイトが役立ちました。
http://hogehoge.tk/tool/number.html
ありがとうございます。
カレー問題
えー、たいへんひさしぶりの更新です。
年度末繁忙期の現実逃避に、結城浩さんのCodeIQのカレー問題をエクセルだけで解いてみました。
解答に書いたメモを転載。本当は実際のエクセルファイルの絵もあるといいのですが、とりあえず。
-
-
- -
-
これならExcel(表計算ソフト)だけでも出来そうだなと取りかかりましたが,30分では終わらず1時間ぐらいかかってしまいました。
- blendlist.txtをExcelで開く。スペース区切りで2列に分け,上に1行挿入して列名(spice1,spice2)を入れる。
- (このままpivot集計をしてみたらspice1とspice2が交換可能であることに気づいた)
- spice1のリスト[A2:A2183]の下のセル[A2184]にspice2[B2:B2183]をコピーして貼り付け,spice2の下[B2184]にspice1[A2:A2183]を追加。
- Excel2003の「データ→ピボットテーブルとピボットグラフ レポート」で新しいワークシート(Sheet1)に行タイトルspice1と列タイトルspice2,項目spice2の個数でピボット集計。対角線に対して対称で,blendlistに存在する組み合わせに1,ない場合は空白の表ができる。spiceの種類は128種であることがわかる。
- ピボットテーブル[A3:DZ133](Sheet1)をコピーして新しいワークシート(Sheet2)の左上セル[A1]に形式を選択して貼り付け(数値)。あとの操作をやりやすくするため。
- 出現頻度が一番高いspiceを探すため,「総計」列の最大値を探すと44が2つある。とりあえず"Minnlabic"を最初の「指標種」とする。
- Sheet1の"Minnlabic"の行(58番目なのでピボットテーブルの行62)を表の下[Sheet1!行135]にコピーして貼り付け。
- spiceの名前[A5:A132]をコピーして[A137]に値貼り付け。
- [B137]に"=B5*B$135"を入力,コピーして[B137:DY264]に貼り付け。[DZ137]に"=sum(B137:DY137)"と入力。コピーして[DZ138:DZ264]に貼り付け。Minnlabicとの組み合わせでポイントになるspiceを共通に持っている数を数える。(28以上と3以下の両極端に分かれるぞ。しめしめ)
- Sheet1の[DZ137:DZ264]をコピーしてSheet2の[EA3]のセルに数値で貼り付け。また[A2]のセルに形式を選択して貼り付け(数値,「行列を入れ替える」にチェック)。
- Sheet2の表[A3:EA130]を選択,EA列の降順で並べ替え。続けて[A1:DY130]を選択。並べ替え→オプション「列単位」で行1の降順で並べ替え。
- この時点で表の対角線に対する"1"の並び方でblendlistが大きく3つのグループになることがわかる。
- "Minnlabic"と相性のよいスパイスを含むブレンドが28以上となるスパイスのグループをグループ1とする。グループ1に含まれないスパイスについて,例えば"Anewatry"を次の「指標種」に選び,同様の処理を繰り返すことが可能。
- しかし,並べ替えた表を見ると,"Minnlabic"と相性のよいスパイスを含むブレンドが0のスパイスの中でLoveaniic, Masema, Peregrama, Rinawatry, Xuomaの5種が,「"Minnlabic"と相性のよいスパイスを含むブレンドが1〜3」のスパイスと相性がよく,残りは別のグループを形成することがわかるので,手作業で5種を"Minnlabic"と相性のよいスパイスを含むブレンドが0のスパイスの中でで最初に来るように並べ替えを行った。(行選択,カットして貼り付け。列選択,カットして貼り付け)
- Group1("Minnlabic"と相性のよいスパイスを含むブレンドが28以上),Group2("Minnlabic"と相性のよいスパイスを含むブレンドが1〜3+上記5種),Group3("Minnlabic"と相性のよいスパイスを含むブレンドが0で上記5種を除く)とすると,それぞれのGroupを含む皿のカレーのブレンドのポイントは以下のようになる。
Group1 | Group2 | Group3 | |
Group1 | 936 | ||
Group2 | 12 | 775 | |
Group3 | 0 | 11 | 448 |
- この際,2皿のポイントの合計を,3つのグループのうちどの1つを別の皿にするかで検討すると,Group3を別の皿にする場合がブレンドリスト中の11の組み合わせだけが達成されないことになり,最大である。
- よってGroup1とGroup2に含まれるスパイスの名前(並べ替えた表の[B2:CQ2])を別のエクセルに貼り付けて昇順に並べ替え,テキストファイルで保存。
- 保存したテキストファイルのタブをスペースに置換。
数値標高モデルデータを実際に使ってみる
オリエンテーリング地図(Omap)は,国際オリエンテーリング連盟の定める規定(http://orienteering.org/resources/mapping/)に従って作成された,ナビゲーションのための地図です。Omapは都市計画図などを基礎図として,トレイル,微地形,林内の通行可能度など,地図を持って山中を走るときに必要な情報を現地踏査により書き加えて作ります。
比較
実際に基盤地図情報の数値標高モデル(http://fgd.gsi.go.jp/download/)を使った地図がOmapの「調査原図」としてどのように使えるか,一つの仮想例としてすでにOmapがあり,各種の数値標高モデル(5mメッシュ標高[航空レーザ測量,写真測量],10mメッシュ標高[1:25000地形図ベース])がそろっている「京都地区」で例示してみます。
Omap「奥大文字・山紫水明東山」 http://www.o-news.net/wr/2010/100530.jpg
この地図を使ったオリエンテーリング大会の記事 http://www.o-news.net/2010/06/z8.php
Omap「奥大文字山紫水明東山」の一部分。磁北方向が上になっているため,下の地図と比較すると約7°左へ回転している。
5mメッシュ標高からGDAL(d:id:kooi:20110713)で5m間隔の等高線を作成し,基盤地図情報縮尺レベル2500の道路縁,建物,水涯線と合わせて表示したもの。等高線が赤茶色の部分が「航空レーザ測量」DEM5A,薄茶色が「写真測量」DEM5Bのデータ部分。基盤地図情報では,同じ地区のDEM5AとDEM5Bは重複しては公開されていない。
基盤地図情報で全国のデータが作成されている10mメッシュ標高から,上記と同様に5m間隔の等高線を作成して表示したもの。
いかがですか。等高線による地形表現をOmapと比較してみると,下図のA,B地点では微小な沢地形が表現されていて航空レーザ測量の精度の高さをうかがい知ることができます。一方,C,D地点ではむしろ10mメッシュ標高の等高線の方が現実の地形に近いと考えられ,写真測量による5mメッシュDEMは(場所によっては)Omapの原図としての精度が10mメッシュ標高より劣ると考えられる場合もあるようです。
数値標高モデルを使った詳細地図
この項はオリエンテーリング地図の作製も念頭に置いて書いています。
基盤地図情報2500には等高線情報は含まれていません。そこで,数値標高モデル(5mメッシュ)から等高線を作成し,基盤地図情報2500とあわせて野外調査の現地野帳として使用できる詳細地図をつくります。ただし基盤地図情報2500はまだ公開されていないエリアも多く,数値標高モデルも10mメッシュは全国で作成されていますが,5mメッシュはまだ一部です。
必要なソフトウェア
GISはオープンソースのQGIS(http://www.qgis.org/wiki/Download)を使用しています(過去エントリ参照)。今回はプラグインのGdalToolsを使用するのでプラグインマネージャで使用可にしておいてください。
数値標高モデルをGeoTiff形式に変換するには,三匹のウリボウさんの「基盤地図情報DEMコンバータVer1.4」(http://space.geocities.jp/bischofia_vb/)を使用させていただきます。
基盤地図情報2500のシェープファイルへの変換は,公共測量ビューア・コンバータ(http://psgsv.gsi.go.jp/koukyou/public/sien/pindex.html)を使います。
データ(基盤地図情報)
基盤地図情報ダウンロードサービス(http://fgd.gsi.go.jp/download/)から,基盤地図情報2500と数値標高モデルの必要な地域のファイルをダウンロードします。
作業手順
1) 基盤地図情報2500は市区町村単位でダウンロードできるので必要な地域をダウンロードし,「公共測量ビューア・コンバータ」でシェープファイルに変換します。
2) 数値標高モデルは(世界測地系)2次メッシュ単位でダウンロードできるので,必要な地域のファイルをダウンロードし,「基盤地図情報DEMコンバータ」で開いて,GeoTiff形式で保存します。
3) QGISで,新しいプロジェクト(座標系はJGD2000)にDEMをコンバートしたGeoTiff形式のファイルを追加し,メニューの「ラスタ」(GdalTools)の「等高線」で,等高線のシェープファイルを作成します。
4) できたシェープファイルをQGISの新しいプロジェクト(座標系は平面直角座標で,オンザフライ座標変換をonにする)に追加します。
5) 必要なシェープファイルを「ベクタ→データマネージメントツール→新しい投影法へエクスポートする」で平面直角座標に変換します。
6) オリエンテーリング地図(OCAD使用 http://www.ocad.com/en/index.htm)へは平面直角座標系に変換したシェープファイルをインポートし,ファイル(レイヤ)ごとに適切な記号をあてます。(測量法で定められた利用申請はしてください)
参考
GDALを使って、基盤地図情報の標高データから等高線のshpファイルを作成する方法 - 自然環境保全のための周辺技術 http://d.hatena.ne.jp/tmizu23/20100215/1266229611
基盤地図情報(標高)をGDALやGRASSで使いたい - 地図はたいへん http://d.hatena.ne.jp/vec2ras/20100107/1262833221
基盤地図情報25000の等高線で計曲線・主曲線・補助曲線を描き分ける
自然環境調査などで使う印刷用3次メッシュマップのQGISを使った作り方(http://d.hatena.ne.jp/kooi/20110619)の補足です。
QGISの地図プロジェクトに読み込んだベクタレイヤの記号は今のところ「共通シンボル」で設定されています。これを「固有値」を使って属性によって記号を変化させるという話です。
基盤地図情報25000の等高線データ(Cntr.shp)の属性ファイル(Cntr.dbf)には「標高」という列があります。この標高フィールドを使って等高線を色分けしてみます。プロパティのシンボルで凡例タイプに「固有値」を選び,分類フィールドを「標高」にします。「分類」のボタンを押すと自動で全ての属性値に個別の凡例をつけることができるようになります。(「標高」のフィールドは数値ではなくテキストでデータが入っています。)
しかしこれを1つ1つ設定するのはかなり面倒です。10,20,30,40,60…と主曲線にあたる標高は「クラスを削除」してデフォルトで一括指定するとしても,50mおきの計曲線だけでも20本以上あります。緩傾斜のところで5mごとに入ってくる補助曲線も何本もあり,それぞれ定義していたらやはり大変です。そこで,ちょうど属性ファイルの「表示区分」フィールドは空欄になっているようなので,ここに等高線の種類を入れて表示区分で凡例を分類することにします。
- QGISの属性テーブルを編集モードにしたときに使える「フィールド計算機」では,まだ複雑な計算は難しいようです。http://www.osgeo.jp/user_guide/user_guidech3.html#x10-1750003.7
- Excel2003まではdbf形式を扱うことができたので,Excelで開いて表示区分のG2セルに以下の式を入れて列へコピーすることが考えられます。しかし,Excel2003までは65535行しか扱えないので,福岡県全部(約29万行)だとオーバーフローします。Excel2007ではdbf形式で保存できないのでこの方法はできません。基盤地図情報閲覧・コンバートソフトでシェープファイルにするときに少しずつ(ダウンロードしたzipファイル1つずつ)別のシェープファイルに変換すればExcel2003でオーバーフローせずにできるでしょう。
=if(right(I2,1)="5"),"2",if(or(right(I2,2)="00",right(I2,2)="50"),"1","0"))
今回は,約29万行のdbfファイルを一発で扱うために,またしても統計解析ソフトウェアRを使います。コードは以下のようになります。dbfはいちおう別名で保存していますので,できたdbfファイルを確認してから,名前をCntr.dbf→Cntr_ori.dbf,Cntr_R.dbf→Cntr.dbfというように変更してください。
library(foreign) CntrType <- function(x) { y <- as.numeric(as.character(x)) ifelse(floor(y/50)==y/50,"1",ifelse(floor(y/10)==y/10,"0","2")) } dd <- read.dbf("C:/GIS/Cntr.dbf") dd$表示区分 <- CntrType(dd$標高) write.dbf(dd, "C:/GIS/Cntr_R.dbf")
Cntr.dbfを置き換えたら,表示区分のフィールドに主曲線は"0",計曲線(50mごと)は"1",補助曲線(5m)は"2"が入っているので,それぞれ適切な記号を定義します。
- R(参考RjpWiki http://www.okada.jp.org/RWiki/)のベクトル演算は強力で,上記のコードではデータの読み込みや書き出しに比べたらあっという間に計算は終わります。関数中でif文ではなくifelseを使ってベクトル演算できるようにしておくのがミソです。
- 個人的には一番手間取ったのはdbfを読み込んだデータフレームで標高の数値がただのテキストではなくてfactorだったので,as.characterをかまさないと期待した結果を返してくれなかったことでした。
インテルメッツォ:地図データを使うのに申請は必要か
国土地理院の地図を使うのには「測量法」できめられた申請が必要ですが,私的な利用等へは申請不要となっています。書籍やインターネットには,限られた量の掲載であれば「出所の明示」で申請不要とされています。
基盤地図情報に関するFAQ。5番目の項目が地図の複製に関することです。
http://www.gsi.go.jp/kiban/faq.html#52
そこでリンクされている「申請フロー」(http://www.gsi.go.jp/LAW/2930-flow.html)を見ると,以下のように書かれています(要約)。
- 私的な使用,教育機関での使用は申請不要
- 学術論文に掲載する場合,刊行物に少量の地図を掲載(挿入)する場合は「出所の明示」で申請不要
- 使用目的により「複製」と「使用」を区別
- 「複製」の場合,サークル内での使用や特定の者へ提出する報告書では申請不要
インターネットへ公開する場合は,300×400ピクセル以下の画像であれば「少量の地図を挿入」と見なされるようです。それ以上1画面範囲内(スクロール不要のサイズ)であれば5枚以内とされています。
http://www.gsi.go.jp/LAW/2930-qa.html
はてなfotofileは画像の長辺が450ピクセルになります。私がこれまであげてきたGISの使い方記事の画像は地図面だけでいえば300×400ピクセル以下のものがほとんどですので,「出所の明示」で申請不要と思っています。