QGISで河川と道路の国土数値情報を重ねあわせる

備忘のためのメモです。参照に来た方、適当になってすみません。

 

1)国土交通省ウェブサイトから、国土数値情報をダウンロードします。

国土数値情報ダウンロードサービス

たいへんお世話になっているサイトです。ありがとうございます!アンケートに答えてダウンロードです。今回は海岸線、道路、河川のデータのうち北海道のものを使いました。

 

2)xmlファイルのものはシェイプファイルに変換。QGISxmlをシェイプに変換するプラグインもありますが、なんだかうまくいかないことも多いです。同じ国土交通省サイトからksjツールをダウンロード(下記サイト)して使うのが確実かと思います。

国土数値情報ダウンロードサービス(JPGIS準拠データ)

3)QGISでベクタ形式ファイルとして順番に読み込む。

今回使用したデータは空間参照系がJGD2000(EPSG:4612)なので、どのファイルもそれで統一する。

4)しかし、河川ファイルは2級河川や準用河川というかなり細かい川までデータになっています。(下記)

f:id:darusmart:20170426211907p:plain

河川データの仕様書を見ると、区間種別コードの指定がありました。属性テーブルで、たとえば1級河川だけに絞れば整理されるということですね。

 

区間種別(COP)コード>

1(1級直轄区間) 2(1級指定区間) 3(2級河川区間) 4(指定区間外) 5(1級直轄区間でかつ湖沼区間を兼ねる場合) 6(1級指定区間でかつ湖沼区間を兼ねる場合) 7(2級河川区間でかつ湖沼区間を兼ねる場合) 8(指定区間外でかつ湖沼区間を兼ねる場合) 0(不明)

 

 

QGISで特定の地物をクリックすると画像が開くようにする方法

QGISオープンソースGISソフトとしてたいへん便利ですが、ときどき思い通りに操作できないことがあります。今回もかなり苦労しました。トライしたのは、特定の地物をクリックして画像を開くという操作です。たとえば、GPS付きのカメラで撮影した画像を地図とリンクさせられたら、とても便利です。またまた苦労を重ね、最後はコンピュータに詳しい同僚の助けをもらって次の方法にたどりつきました。

1)画像からGPS情報を抽出する。

まず、GPS情報を持っている写真の情報(EXIF)を一括して閲覧したり、エクスポートしたいと思いました。枚数が少なければ手入力でもよいですが面倒ですね。これにはExif情報を抽出するフリーウェアを利用しました。リンクを掲載してよいかどうか不明でしたのでここでは名前は伏せます。私のカメラ(Canon PowershotD30)では緯度経度が「・・度・・分・・秒」となっていたのでこれをExcel上のワークシート関数で変換しました。たとえばC2に上記の緯度がある場合、次の関数で小数点に変換できました。

=LEFT(C2,FIND("度",C2)-1)+(MID(C2,FIND("度",C2)+1,FIND("分",C2)-FIND("度",C2)-1))/60+(MID(C2,FIND("分",C2)+1,LEN(C2)-FIND("分",C2)-2))/3600

今回はヨーロッパの写真が対象で、東経と西経の両方が登場しました。区別はそれぞれプラス、マイナスとしておけば大丈夫でした。これをCSVファイルにして、QGIS(今回は2.12.0を使用。PCはwindows7です)から読み込みます。いつもどおりセルやファイル名に日本語を含めないのが無難です。このCSVファイルに画像の場所を指定する列も加えておきます(列の一部はこちらの都合で隠しています)。

f:id:darusmart:20151120174123p:plain

 

QGISからのCSVファイル読み込みは「レイヤ>ベクタレイヤの追加>デリミティッドテキストレイヤの追加」です。

2)QGISのアクションとバッチファイルを利用して地物と画像のリンクをつくる。

QGISにはプロパティ(レイヤーを選択してダブルクリックすると開くメニューの中にある)を開くと、”アクション”というタブがあります。このアクションを利用すると地物と特定の動作を結びつけることができます。ところが少しクセがあるようで、いくつかのウェブサイトに書かれていた方法ではうまくいきませんでした。同じように困った方は次を試してください。

(1)バッチファイルをつくる。

次のようにテキスト入力したファイルに名前をつけて、拡張子をbatとして保存します。フォルダはプロジェクトと同じところに入れておくのが無難です。

f:id:darusmart:20151120185118p:plain

(2)アクションを次のように設定する。

下のほうにある欄で設定をします。”タイプ”は一般、”名前”は表示のためだけなので日本語でも大丈夫です。アクションのところにはさきほど指定したバッチファイルの場所を入力します。そのあとスペースをはさんで画像の場所を記入したカラムを示します。アクションの欄の下にあるポップアップで選んでから”フィールドを挿入”とします。最後に忘れず”アクションリストへ追加”をクリックします。

 

f:id:darusmart:20151120190456p:plain

上記によってアクションリストに下記が加わりました。OKを押して通常の画面に戻ると、下図で地物アクションの実行というアイコンのポップアップに上記アクションが追加されるはずです。ここでアクションを選んだら、あとはそれぞれの地物をクリックすると画像が開くはずです。(OSMをベースマップにしています。)

f:id:darusmart:20151120191357p:plain

バッチファイルを使わなくてもできるようですが、私は今のところうまくいっていません。画像を開くとき、一瞬バッチファイルが開くのがうっとおしいかもしれませんが、とりあえず目標は達成しました。

 

 

 

 

 

 

 

Rでヒストグラムをつくる

ヒストグラムは統計の授業でも最初のほうに出てくる基本的なグラフですが、Excelでつくろうとするとアドインを加えたり、区間を定義したりとちょっと面倒です。Rだとかなり簡便に作れます。

1)データの読み込み

csvファイルを読み込ませる方法もありますが、ウェブで調べてみるとコピペでデータを読み込ませる方法があり、これがたぶん一番簡単ではないかと思いました。まず、ヒストグラムをつくりたいデータをコピーします。Excel上に入力してあるときはデータ範囲をドラッグしてコピーです。次に、Rのコンソール画面で次のように入力します。

>A=scan("clipboard")

これでEnterを押すと、さきほどのデータが読み込まれます。たとえば140個のデータなら、

Read 140 items

などと返されるでしょう。Aというのはデータセットに名前をつける便宜的なものなので"D"でも"S12"でも"Money"でもなんでもよいです。

2)ヒストグラムの作成

> hist(A)

とするともうヒストグラムが表示されるはずです。

3)区間の変更

デフォルトで設定されている区間でなく指定したい場合は次のように調整します。

> hist(A,breaks=seq(0,10,0.5))

これは0から10までの範囲で0.5ずつの区切りでデータを分けたヒストグラムを作成せよという意味になります。

4)データの重ね合わせ

別のデータを上記ヒストグラムに重ねるときは、まず別のデータセットを読み込んでおきます。上記1と同様にコピペです。

>B=scan("clipboard")

 つぎにBのヒストグラムを重ねわせるには、

> hist(A,breaks=seq(0,10,0.5))

> hist(B,breaks=seq(0,10,0.5),col="#808080",add=T)

とします。col="...."は色指定をするコマンドで、色のコードは6ケタで与えます。このコードは次のサイトを参考にしてください。add=Tが先にできたヒストグラムに新たなヒストグラムを重ね合わせるコマンドです。

カラーコードで色指定

y軸がおさまらないときは、上記コマンドの中にylim=c(0,100)等と入れます。y軸の範囲が0から100までという意味です。たとえば次のように書きます。

> hist(A,breaks=seq(0,10,0.5),ylim=c(0,100))

 

CONTAXレンズを使ってみました

フィルムカメラの時代はもうすっかり過ぎ去った感があります。私は10年ほど前にフィルムカメラNikon FM3Aを購入しましたが、今となってはほとんど使用していません。代わってもっぱら使用しているのはOLYMPUSデジタル一眼レフカメラE-3です。耐久性も高くて機能も申し分ない良いカメラです。しかし、レンズには不満がありました。もちろんOLYMPUSの純正レンズにも良いものはあるのですが、デジタル対応はけっこう良い値段がします。一方、フィルムカメラ用のレンズは中古市場でどんどん値下がりしてきました。そんなわけで、今回中古でCONTAXのTele-Tessar(テレテッサー) T*300mm F4 MM というレンズを購入しました。光学系はツァイス社製で、なんと発売時は21万円の値がついていたこのレンズが、今回3万円台で購入できました!もちろんアダプターをかまさなければOLYMPUS E-3のフォーサーズマウントには使えませんし、ピントも露出もマニュアルです。きっとデジタル世代には考えられないことでしょうが、30年以上前からフィルムカメラを使っていた身には懐かしい感覚です。もちろん動き回る被写体は不利ですが、あわせたいところにピントをもってくるとか、露出を自身の希望でこきざみに変えるにはかえって楽なくらいです。しかもカメラはデジタルなので結果は確認できる(笑)。かつては現像結果をドキドキしながら見たものですが、そのドキドキ感は良い思い出として、デジカメで確認できることによって実はマニュアルレンズはけっこう使えます。撮影した4つの写真を添付しました。シャープさと上品さを兼ね備えている印象です。1kg近くとそれなりに重いですが、シャッタースピードに気をつければ手持ちは十分可能です。

f:id:darusmart:20150507204708j:plain

f:id:darusmart:20150507204714j:plain

f:id:darusmart:20150507204729j:plain

f:id:darusmart:20150507204724j:plain

名古屋でおすすめの博物館

家族旅行で名古屋、岐阜にでかけました。名古屋でトヨタの博物館を観ることにしました。さすがに世界のトヨタです。広大な面積の施設に、自動車や紡績の工場がそっくり入ったような博物館がありました。綿がよじられて繊維になってゆく仕組みを実際に体験しながら学べます。子供には携帯ストラップや輪ゴムで走る車の模型プレゼントもありました。これで大人1人500円は安すぎです!

www.tcmit.org

f:id:darusmart:20150330004317j:plain

f:id:darusmart:20150330004400j:plain

弘前市の果樹園と水田の分布図を作成する

業務で、弘前市のりんご畑と水田の分布図を作成することとなりました。

電子情報の整備された今なら簡単にりんご畑の図くらい書けるだろうと甘くみていました。リンゴ畑のデータがなかったのです。ならば、果樹園ならすぐに手に入るだろうとおもいました。こちらも市内での分布を示すような細かいメッシュデータがありません。そこで、ふと思いつきました。「2万5000分の1地図に出てくるような果樹園マークをポイントにしていけばいいのでは?」ということです。弘前市の果樹園はりんごが主体であることは間違いありません。ほかにも含まれるものの、今考えられる方法はそれしかありませんでした。

次は作業手順です。まず、QGIS上に基盤地図情報から弘前市を表示させました(これは国土地理院のウェブサイトやQGISの解説本を参照ください)に説明があります。そして、ひたすら果樹園マークをポイントデータにしてゆきました。

f:id:darusmart:20150222165121p:plain

ポイントデータはシェイプファイルの1種です。QGISメニューのレイヤー>レイヤーの作成>新しいシェイプファイルレイヤーと選び、点(ポイント)が選ばれていることを確認します。あとは、点を打ってはidを順番に入れてゆくだけです。100,200,300・・・と、このあたりで自身の計画の無謀さに気付き始めます。「けっこう多いな・・・」。「1000、2000・・・」これは危険な作業だと気付いたときにはすでに時遅しです。途中で保存し忘れて400ポイント失う事件などもあり、その後は100ポイントずつで保存するようにしました。のべ1週間の作業となり、最終的12,000ポイントとなりました。弘前には一面のリンゴ畑が広がっているのです!

そして、水田は既存データから入手できます。国土交通省のウェブ(下記)で提供されている土地利用細分メッシュデータです。

国土数値情報ダウンロードサービス

しかし、このデータはxml形式でそのままだとQGISでは読めません。そこでxmlからshpへの変換ツール(ksjツール)も提供されていたのであわせてダウンロードしてそこから変換を試みました。・・・しかし、うまくいきませんでした。QGISプラグインで”fgddemImporter”というものが”Minoru Akagi”さんによって提供されています。

fgddemImporter プラグイン

こちらを利用したところうまく読み込めました。

しかし、土地利用細分メッシュデータは青森県で4つくらいと大きなデータです。そこで、弘前市の範囲を指定する必要があります。まず、先ほどと同じ国土数値情報のサイトから、行政区域を選択して青森のデータをダウンロードします。そして、やはりfgdemImporterで・・・と思ったら、今度はうまくいきません。こちらはksjツールで変換できました。なぜだか私にはわかりません。

土地利用細分メッシュデータも弘前市以外の部分をまず削除します。レイヤを選択してから、ベクタの編集モード(えんぴつマークをクリック)にして、範囲を選択します。そして、選択物の削除を選ぶと削除できます。削除する量が多いとフリーズしたように見えますし、時としてほんとうにフリーズします。適度な量ずる削除しましょう。

f:id:darusmart:20150222172533p:plain

 

f:id:darusmart:20150222172753p:plain

 

 

市町村データのほうも同様です。読み込んだ青森県各市町村のデータのうち弘前市だけを抜き出します。シェイプファイルなので、属性テーブルを開いて弘前市以外のセルを削除してゆけばよいです。

 そして、弘前市の行政境界で水田のデータを切り抜きます。これは

ベクタ>空間演算ツール>クリップと選び、

まず水田のデータ(入力ベクタレイヤー)を選びます。つぎにクリップレイヤーとして弘前市の行政データ(ポリゴン)を選びます。すると、弘前市のみの水田データになります。

完成した弘前市のリンゴ畑(ピンク色)と水田(黄緑)マップです。なにも複雑なことはしていませんが、実にたいへんな作業でした。

  f:id:darusmart:20150222173448p:plain

 

 

 

Rでコルモゴロフ・スミルノフ検定

やや古典的なひびきの統計手法ですが、分布の正規性や、2標本の代表値や分布の違いを検定するときに使います。2標本の場合は代表値(一般には平均値でしょうか)が同じなら分布の検定になるという関係だそうです。詳細はほかのウェブにもっと正確な文章が出ているので、ここではRでの検定手法だけメモしておきます。

1)1試料のとき、ここでは正規性の検定をする場合を書きます。

まず、データの読み込みです。「ksdata」という名前でこれから検定するデータを読み込ませます。Rを起動して直接入力する場合です。

>ksdata<-c(90,126,232.41,94.5,156,120.9,180.4,32.25,29.2,100.11,21,23.1,77,78.75,122.32,59.78,95.59,124.95,112.8,98.4,31.39,31.05,26.4,184,72.8,62.98,23.56,53.35,49.5,33.48,99.33,39.2,11.34,31.92,41.8)
このようにカンマ区切りで指定すればよいです。

続いて検定。

> ks.test(ksdata,"pnorm",mean=mean(ksdata),sd=sd(ksdata))


pnormが正規分布のときに使用するコードです。正規分布は平均と

標準偏差(sd)を指定すれば形が書けるので、それぞれksdataのものを使うのだと指定しています。結果は次の通り。p値が0.567なので、有意水準5%とすれば正規分布といってよいかと思います。


        One-sample Kolmogorov-Smirnov test
data:  ksdata
D = 0.1284, p-value = 0.567

 

2)2試料のとき

上記のksdataと比較するデータセットをもうひとつ読み込ませます。ks2dataと

しました。

>ks2data=c(94.4,117,230.4,91,166.32,139.12,206.48,56.76,29.6,96.75,19.2,23.18,104,82.5,115.2,65.52,100.33,113.9,132.8,96,29.4,49.92,26.8,149.73,74.52,61.36,23.4,39.22,49.02,49.02,33,102,33.44,34.04,34.79)

2試料を比較します。

> ks.test(ksdata,ks2data)

結果は次の通り。

 

 Two-sample Kolmogorov-Smirnov test

data:  ksdata and ks2data
D = 0.1143, p-value = 0.9763
alternative hypothesis: two-sided

 警告メッセージ:
In ks.test(ksdata, ks2data) :
   タイがあるため、正確な p 値を計算することができません

タイとは同一値のデータとのことです。なんだか不安を感じるメッセージがあるものの、とりあえずp=0.9763なので2つの分布に差はないとしてよいかと思います。それぞれの平均も下記のとおりです。実はクロテンの左右の精巣サイズのデータでした。片方が人間だと少し下がっているとかいいますが、たぶんサイズはそれほど違わないのではないでしょうか!

> mean(ksdata)
[1] 79.07314
> mean(ks2data)
[1] 82.00343