QGIS:OpenLayers Pluginの背景画像をコンポーザで印刷するとずれる

タイトルの質問がとどいたのでちょっとしらべてみました。
2年前から問題が指摘されていたみたい。
http://hub.qgis.org/issues/5827
#13のコメントの内容を再現できたのでメモしておきます。

コンポーザでたとえばA4サイズで以下のような地図を印刷するとします。
以下の画像のように余白を設定することはわりとよくありますね。
p001
そうすると重ねてあるベクタレイヤがずれて印刷されてしまいます。
質問の状況が再現できました。
p002
これはアイテムプロパティの”位置とサイズ”が以下のように
印刷サイズであるA4のサイズと食い違っているために発生するようです。
p003
印刷地図の余白がなくなってしまいますがとりあえず
“位置とサイズ”をA4サイズと同じ値に修正し印刷してみます。
p004
ずれは修正されているようです。
p003
OSMの地図を背景にした印刷需要は高そうなので
ちょっと不便だなーと思っているヒトもかなりいるだろうな。

航空レーザ測量データをひたすら ESRI shapefile に変換するだけ

ディレクトリに格納されている複数の航空レーザ測量のデータ(xxxx.txt)を
さっくりESRI shapefile形式(ポイント)に変換したいという相談があったので
pythonで書いてみました。

カンマで区切られたテキストの2列目がx 3列目がy 2列目がz
となっていたのでそれにあわせて書いてあります。

ogr2ogr使えばいいんだけど
筋トレ的に書いたので気にしない気にしない。

pyshpというライブラリが必要です。
コマンドプロンプトから
pip install pyshp
または
easy_install pyshp
でインストールできます。

pip も easy_install もインストールされていない場合は
このあたりからバイナリ(pyshp-1.2.0.win-amd64-py2.7.exe)をダウンロードして
インストールしてください。

その後以下をテキストエディタ等で適宜名前(ex. lp2shp.py)をつけて
Python27\Lib
直下に保存してください。

# 10/9追記:余計な処理をしてメモリを無駄に消費していたので修正



その後、python(command line) 等で



とテキストデータの入っているディレクトリを指定して実行してください。
同じディレクトリ内に xxxx_p.shp というファイル名でshapefileが出力されます。

カメラトラップの撮影記録を撮影間隔&種でグループ化

前回の記事の続きです。

以下のような自動撮影装置の撮影記録 CSVファイルがあって



X番目に撮影された動画とX+1番目に撮影された動画の撮影時刻を比較し、
差が指定された時間以内(今回は分単位)
かつ
同種同性の個体が撮影されていれば
両者をグループ化する(末尾の列に同じID(gid)を追加する)というスクリプトをpythonで書いてみました。
データ処理後以下のような出力結果が得られます。



以下の通りメモとして残しておきます。



QGIS地図帳機能テンプレート(シンプルな帳票作成)

QGISの地図帳機能を使って
帳票を出力する作業をアルバイトさんにお願いするために
シンプルなテンプレートデータを作りました。

ベクタデータの属性テーブルから各種の情報を呼び出して
プリントコンポーザ図面内の要素として出力する、という手順を
テンプレートをいじって覚えていただこうということです。

いまのところアルバイトさん2人しか使わないのも
なんだかもったいないので興味のある方はここからダウンロードして
お試しくださいませ。

とりあえず動かすには以下の手順で操作してみてください。

なおこのテンプレートですが Windows かつ QGIS2.4 でのみ動きます。

(MacだとHTMLの記述を行うところが機能しません…
詳しい方なぜだか教えていただけるとうれしいです)

1) 解凍した atlas_sample フォルダをC:\temp\に保存する
(“そんなとこイヤ”というかたはとりあえずどこかに保存して進めてくださいね)

2) QGISを起動しタイルレイヤプラグインをインストール
プラグイン→プラグインの管理とインストール
→TileLayer Pluginを選択→プラグインのインストール

3) タイルレイヤプラグインの設定を行う
web → タイルレイヤプラグイン → タイルレイヤを追加する
→ 設定 → 外部レイヤ定義ディレクトリ
→ 同梱している “tile_info” フォルダを指定 → OK

4) atlas_sample.qgs を開く

5) プリントコンポーザを開く
プロジェクト→コンポーザマネージャ→atlas_sampleを選択→表示

++++”そんなとこイヤ”なひとはここを変えてね+++++

現場写真と書かれている場所の下にある
“ここをクリック”をクリック→右側にある”アイテムプロパティタブ”を開く
ラベルのメインプロパティ内4行目の
C:/temp/atlas_sample/pic/とあるところを
保存した場所(ex. D:/hoge/atlas_sample/pic/)へと書き換えてください。

++++++++++++++++++ここまで++++++++++++++++++++

6) 地図帳の表示と出力
地図帳→地図帳のプレビュー で出力される地図が表示されます (まだ写真は出ません)
次の地物へをクリックすると場所が変わります。

7) 地図帳をPDFへと出力する
地図帳→地図帳をPDFへと出力する
→出力場所とファイル名を指定→保存 (ここでやっと写真が表示されます)

なお 背景図として 国土地理院 地理院タイル 標準地図 を利用しています。

注意:写真はあくまでサンプルデータであるため現場の状況とは一致していませんよ!

ざっくり内容を説明すると…

++データの構造++

・仮想の調査地点を示すベクタデータ(ポイント)(sample_location.shp)
・調査地点番号と対応したCSVテーブル(loc.csv)

の2つのデータのみです。

「sample_location.shp」に「loc.csv」を
地点番号をキーに結合させて使用しています。

++地図帳を使うには++

プリントコンポーザ画面の右ペインにある「地図帳の作成」 タブを開いて

・地図帳の作成にチェックを入れる
・被覆レイヤに調査地点レイヤ(ここでは sample location)を指定する
→地図帳ツールではこれを中心に配置するように表示が切り替わります。

++属性フィールドの情報を図面に挿入するには++

・ラベルを追加してプロパティの操作を行う
→緯度経度など図上の各アイテムを触って「アイテムプロパティ」
のラベルのメインプロパティを表示させてみてください。

直接そこに書き込む or 式の挿入を利用して
画面表示要素を作成できます。

→写真の挿入については
++”そんなとこイヤ”なひとはここを変えてね++を参照してください

かなりテキトーな説明ですみません。
ではでは。

カメラトラップ(camera trap)の画像データを撮影間隔でグループ化

地理情報とあまり関係ないですが
先日参加した学会の懇親会にて主催者の一人が
“こんな処理ができたらいいなー”
とおっしゃっていたことを次の週の野外調査の宿舎で思い出しました。
んでビールのみながらPythonでスクリプトを書きました。

自動撮影装置(自動撮影カメラ)で
X番目に撮影された画像(動画でもよい)とX+1番目に撮影された画像の
タイムスタンプを比較し、差が指定された時間以内(今回は分単位)であれば
両者をグループ化する(同じフォルダに格納する)という単純な処理です。



とりあえず
ここからダウンロードして解凍
gpic フォルダをフォルダごと
windowsであれば
C:\Python27\Lib\
にコピーしてください。

それで

python (command line)
または
コマンドプロンプト(最初にpythonと入力しENTERする)

などで以下のようにコマンドを入力して実行してみてください。



対象とするデータフォルダ内にグループ毎のフォルダが出力され
そこに対象グループに含まれるファイルが格納(コピー)されます。

取り急ぎメモかわりに書きました。
ではでは

TOPEXをGRASSで算出してみる

もくもく単純データ処理に嫌気がさしてきたので
気分転換を兼ねて,
以前相談を受けていて放置していた
TOPEX(Topographic Exposure)の算出を
Grassをつかってやってみました。

GrassでTOPEXを計算,ていうのはこの方がやっていたので参考にしました。
Calculating Togographic Exposure With Grass

もともと俯角はすべて0°とするものだったらしいけど
負の値として計算することもあるようですね。
たとえばこれ(P47-48)
EMPIRICAL MODELLING OF WINDTHROW RISK USING GEOGRAPHIC INFORMATION SYSTEMS

DEMを基にしたTOPEXは
こんな感じで計算するようです。

1) 対象セルA00から8方位の計測線を延ばす
2) その延長上かつ特定の計測範囲内において
一定の間隔でA00からの高低角を算出する
(下図では A00-A33 , A00-A55
およびA00-A77 間の高低角を算出する)
3) 求めた高低角の最大値を
その方位の高低角と定義する
4) 各方位で同様の計算を行い
8方位すべての値を積算したものをTOPEXとする

↓かなりアレな図ですが(平面図です)背景のメッシュはDEMだと思ってください…

tekitou

どうやら地上開度とアイディアはほとんど同じ。
開度による地形特徴の表示

とりあえず

1)Grass を外部からうごかす環境設定済み
2)DEMはGrassにインポート済み
3)俯角もそのまま計算する(負の値として処理)

という前提で進みます。

引数多すぎとか8方位なのにいらん計算してるとか
そもそもGeotiffのインポートからやれよとか
いろいろあるけど気分転換なので 気にしない 気にしない

ここからダウンロードして解凍
topog フォルダをフォルダごと
windowsであれば
C:\Python27\Lib\
にコピーしてください。

それで

python (command line)
または
コマンドプロンプト(最初にpythonと入力しENTERする)

などで以下のようにコマンドを入力して実行してみてください。



topex.pyの中身はこうなっています。

環境変数の設定とかめんどくせ という場合は
同梱してあるtopexf.pyを使って

1)出力コマンドをテキストファイルに出力
2)grassのレイヤマネージャにあるコマンドコンソールにコピペして実行

てのがてっとりばやくて楽かもしれないですね。
これはpython2xがあれば動きます。



topexf.pyの中身はこうなっています。



コメントはかなり適当(ごめんなさい)
これみたいにGDALでやったほうがべんりなんだろうなー
gdalでfocal statistics

ではでは。

属性で ベクタデータ をサブセットする

昨日夜10時
「明日までにxxフォルダ内の等高線ベクタデータから100m間隔のデータを作っといて」
とさっくり依頼された。
数百ファイルにわかれて30GBくらいあるんだけど・・・

finds.jpの WMSじゃダメ? と一応聞いてみたが
ローカルに持っておきたいとのこと。

しょうがないので以下の内容(作者さまに感謝)をちょっと変えて
3の倍数でアホになるアレで以下のような処理を行った。
https://code.google.com/p/geospatialpython/source/browse/trunk/subset.py

あ、でも OGR を使えばもっと簡単だったのかも。



外部からGRASSの機能を呼び出したいの

++ 以下よりFOSS4Gアドベントカレンダー2013への投稿記事となります ++

最初からタイトル間違ったんじゃないかと疑いたくなる話題ですが

  • >>FOSS4G 勉強会 in 函館

11/30-12/1 と2日間にわたり
FOSS4G 函館勉強会
(於 北海道大学函館キャンパス)
に参加してきました。
まずは会場の雰囲気をお届けしますね。

IMG_2298_coreday

1日目の事例紹介セッションの様子です。
質問がよく出ていました。

IMG_2312_handson

これは 2日目 QGISハンズオン
の様子です。
主要な参加者は学生さんだったのですが
地元企業の方もいらっしゃっており
道南ユーザのよい顔合わせとなったのではと思います。

さて私は1日目の事例紹介セッションで
お話させて頂きました。
内容をどうしようかと考えた挙句
みんな QGIS QGIS ってチヤホヤしやがって・・・
との嫉妬心からGRASSの話にしました(ウソ)

  • >>外部からGRASSの機能をつかおう

客層と内容が合致してなかった残念なプレゼンはさておき(さておくなよ)
発表の中ででもちょっとだけ紹介しました
外部からPythonを介してGRASSの機能を使用する
ことについて書きます。
ここでやっとタイトルの内容です。

エクストリームなpython使いの方には
いまさらな話だとおもいますが
取り組みだした当初、日本語で解説したサイトが見つからず
ちょっと苦労したこともあり、気にせず書くことにします。

以降

OSがWindows
OSGeo4W経由でGRASS6.4.3をインストール済み

ということを前提に進めていきますね。

内容まちがってる!てのにお気づきのかたは

ぜひぜひご指摘いただければと思います。

  • >>システム環境変数

最初にシステム環境変数を追加・編集する必要があるので
それについて。

Pythonの環境変数PATHの設定がまだの場合は
以下の場所の情報を参考にして行ってみてください。
http://www.pythonweb.jp/install/setup/index1.html

なお上記の設定中に出てくるフォルダ
Python26はPython27と読み替えて設定願います。
インストールしているソフトウェアによっては
より下位のディレクトリにpython.exe
が入っていることもあります。
それに伴って設定するpathも変更してください。

さてシステム環境変数を
新規追加または既存のものを編集していきます。

以下より

変数名 変数値

の順番で書いていきます。
なお単一の変数名に対して変数値が複数にわたる場合は
セミコロンで区切って入力してください。



こんなところです。
ちょっと多くて面倒ですね。
(マシンを換えるたびに あれ、どーやるんだっけ と思いメモを検索してます)

  • >>スクリプトを書くよ

実際にスクリプトを書く際には



てなのを最初に書いておき
あとは自分の行いたい処理を書いていきます。
自分のあまりに具体的な処理を載せるのもアレですので
あっさりとした繰り返し処理のサンプルを書いておきます。



これでスケール間比較を行いたいから
複数(たくさんの)パターンのデータを作って欲しい
といういかにもありがちなリクエストにも
お答えできるやもしれません。

電子国土基本図のオルソ画像をQGIS上で表示する

元ボスと打ち合わせの最中に、
『電子国土オルソ画像をQGISで表示できないかな?』
と聞かれた。
以下の記事が頭にあったので

地図はたいへん QGISで電子国土基本図

『地図画像ができるので オルソもできますよ。たぶん。』
とお答えしたら間髪いれずに
『やって』
・・・ハイ。ではやりますね。

完全にひとのふんどしでなんとやらな内容なので気がすすまないのですが
検索してもオルソの表示方法に関しての具体的な記載がみつからなかったので
私をふくめた一般ユーザ向けとしては意味はあるのかなあと思います。

以降

地図はたいへん QGISで電子国土基本図

に記載されている設定をすでに行っている前提で進みます。

まず電子国土ポータルの
電子国土Webシステム(固定ピクセル版)(Ver.4)OpenLayersサンプル集
をのぞいてみる。

そこで デフォルトで地図表示を行うために必要な記述は



という2行らしいことがわかった。

次にオルソを表示するために必要な記述をみてみる。



オルソを出すには地図表示のための2行に加えて
・空中写真データセットを取得する getOrthoDataSet :2行目
・取得したものをデータセットとして設定する setDataSet :3行目
という2行を足せばよいらしい。

そして

地図はたいへん QGISで電子国土基本図

の記事にもどって cjp.html をみてみる。



地理院のサンプルと同じ記述があった。
ではここに2行加えればよいのかな。

4行のつじつまはあっている。
これでたぶんオルソを表示させられるはず。

この2行を追加したたファイルを cjp_ortho.html とし
($HOME)/.qgis/python/plugins/openlayers/html に保存しておく。

+++++++++++++++++
最近openlayers pluginが1.1.0にアップデートされました。
パスがほんの少しだけ変わっています。
($HOME)/.qgis/python/plugins/openlayers_plugin/html

QGIS2.0にアップグレードされたかたは以下の場所になります。

($HOME)/.qgis2/python/plugins/openlayers_plugin/html
+++++++++++++++++

cjp_ortho.htmlの中身は以下の通り



空中写真のデフォルトデータセットでは電子国土基本図(オルソ画像)と呼ばれる2007年以降に撮影された画像を指定している。これらの画像の撮影範囲は平野部かつ都市域に限られている。

そこで全国網羅されている国土画像情報第1期(1974年~1978年)版も別につくっておく。

以下の地理院サイトの記述から、ズームレベル毎に背景地図データを指定すればよいこと、また、背景地図はdataIDというもので指定することがわかった。

電子国土Webシステム(固定ピクセル版)(Ver.4)について
背景地図タイルの呼び出し方法

電子国土Webシステム(固定ピクセル版)(Ver.4)試験公開時の背景地図データ
Ver.4(試験公開)の配信背景地図一覧

1970年代の空中写真オルソのIDは NLII1 である。

上記をふまえてcjp.htmlを以下のように変更し
cjp_ortho70s.html としてさきほどと同じディレクトリに保存しておく。



var dataSet 以下のdataIdのうち, ズームレベル15~17をNLII1に指定し、
その他はデフォルトデータセットと同じとした。

次に openlayers_plugin.pyの編集に移る。

($HOME)/.qgis/python/plugins/openlayers にある openlayers_plugin.py をひらく。

+++++++++++++++++
↑ここも同様にパスが変わっています。
($HOME)/.qgis/python/plugins/openlayers_plugin

QGIS2.0にアップグレードされたかたは以下の場所になります。

($HOME)/.qgis2/python/plugins/openlayers_plugin

+++++++++++++++++

# Layers のブロックの下に以下の2行を挿入



これを保存し同じディレクトリにある openlayers_plugin.pyc を削除しておく。

QGISを起動し メニューバーの プラグイン→ OpenLayers plugin をひらくと
電子国土のオルソ画像への選択肢が追加されている。
ためしに開いてみたところ、うまく表示されているようだった。

2012年の始まり

1月10日の朝、出勤すると机の上にこれが置かれていた。

年始の挨拶が書かれた付箋が付いていたのでどうやらプレゼントらしい。

このメモの存在は知っていたが実際に手にしたのははじめてだ。
何らかの理由で使用されなくなった地形図を切断し裏向きに綴じたもので、
メモ面の裏にいろいろな地域の地形図の一部をみることができる。

最初のページをめくってみると・・・

北海道内の図幅の一部だ!(増毛町)

これがきっかけで、メモ帳にはいったい何枚北海道内の図幅が含まれているかを知りたくなり
80ページすべての裏をチェックしてみた。

いうまでもないが、お仕事という現実からの逃避行動である。

北海道内/北海道外の判別には河川名や市町村名を用いた。
また、”カリベツ”や”ペーペナイ”などのアイヌ語名称が入っているのも
重要な判別ポイントである。

結果、80ページ中18枚が北海道内の図幅に含まれるものであった。

へー結構多いなあ、18/80だから22.5%かあ、とか適当に暗算していたが
ふと思い立ち、北海道の総面積と日本の総面積を調べてみた。
それぐらい知っとけ自分・・

北海道:83,457km2
日本:377,914km2

(あくまで地形図が刊行されている地域の面積である)

実際の面積比とおなじくらいだ!

・・・そんなどーでもよい発見とともにはじまった2012年、

皆さま今年もよろしくお願いいたします。