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

投稿日時: 2014年11月8日投稿者: admin

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

コンポーザでたとえばA4サイズで以下のような地図を印刷するとします。
以下の画像のように余白を設定することはわりとよくありますね。

そうすると重ねてあるベクタレイヤがずれて印刷されてしまいます。
質問の状況が再現できました。

これはアイテムプロパティの”位置とサイズ”が以下のように
印刷サイズであるA4のサイズと食い違っているために発生するようです。

印刷地図の余白がなくなってしまいますがとりあえず
“位置とサイズ”をA4サイズと同じ値に修正し印刷してみます。

ずれは修正されているようです。

OSMの地図を背景にした印刷需要は高そうなので
ちょっと不便だなーと思っているヒトもかなりいるだろうな。

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

投稿日時: 2014年10月7日投稿者: admin

返信

ディレクトリに格納されている複数の航空レーザ測量のデータ(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追記：余計な処理をしてメモリを無駄に消費していたので修正

import shapefile
import csv
import glob
import os.path
def lp2shp(fd):
    # ファイルリストの作成
    flist = glob.glob(fd + "*.txt")
    # ファイルリスト毎に繰り返し処理
    for i in flist:
        # テキストファイルの読み込み
        cobj = csv.reader(open(i, "rb"))
        # ESRI shapefile(point)をつくる
        w = shapefile.Writer(shapefile.POINT)
        # 属性フィールド elev(標高値)を float, 出力フィールド幅:15, 精度:5, で作成する
        w.field("elev", "F", 15, 5)
        # 属性フィールド flagを integer 出力フィールド幅:5 で作成する
        w.field("flag", "N", 5)
        # テキストの行毎の情報をもとにポイント作成        
        for j in cobj:
            # 2列目をx座標 3列目をy座標としてポイント作成
            w.point(float(j[1]),float(j[2]))
            # 最初の属性フィールド(elev) に テキストファイル 4列目の情報を入力
            # 2番目の属性フィールド(flag) に テキストファイル 5列目の情報を入力
            w.record(float(j[3]),int(j[4]))
        # テキストファイルのフルパス情報からファイル名のみ抽出
        os.path.basename(i)
        # ESRI Shapefileの保存(テキストファイルのファイル名に"_p"を追加した名称を付与)
        w.save(str(i)[0:-4] + "_p.shp")
    print "finished"

if __name__ == '__main__':
    print "this is code block"

import shapefile

import csv

import glob

import os.path

def lp2shp(fd):

# ファイルリストの作成

flist = glob.glob(fd + "*.txt")

# ファイルリスト毎に繰り返し処理

for i in flist:

# テキストファイルの読み込み

cobj = csv.reader(open(i, "rb"))

# ESRI shapefile(point)をつくる

w = shapefile.Writer(shapefile.POINT)

# 属性フィールド elev(標高値)を float, 出力フィールド幅:15, 精度:5, で作成する

w.field("elev", "F", 15, 5)

# 属性フィールド flagを integer 出力フィールド幅:5 で作成する

w.field("flag", "N", 5)

# テキストの行毎の情報をもとにポイント作成

for j in cobj:

# 2列目をx座標 3列目をy座標としてポイント作成

w.point(float(j[1]),float(j[2]))

# 最初の属性フィールド(elev) にテキストファイル 4列目の情報を入力

# 2番目の属性フィールド(flag) にテキストファイル 5列目の情報を入力

w.record(float(j[3]),int(j[4]))

# テキストファイルのフルパス情報からファイル名のみ抽出

os.path.basename(i)

# ESRI Shapefileの保存(テキストファイルのファイル名に"_p"を追加した名称を付与)

w.save(str(i)[0:-4] + "_p.shp")

print "finished"

if __name__ == '__main__':

print "this is code block"

その後、python(command　line)　等で

import lp2shp as l
l.lp2shp("C:\\temp\\")

1 2	import lp2shp as l l.lp2shp("C:\\temp\\")

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

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

投稿日時: 2014年10月1日投稿者: admin

返信

前回の記事の続きです。

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

ID_A,地点名,ファイル名,日,時,日時,種名,ID_B,雌雄,行動
1,A,A1.avi,2014/01/01,10:00,2014/01/01 10:00,エゾシカ,3,メス,採餌中
2,A,A2.avi,2014/01/01,10:05,2014/01/01 10:05,エゾシカ,3,メス,採餌中
3,A,A3.avi,2014/01/01,11:00,2014/01/01 11:00,エゾタヌキ,3,オス,採餌中

ID_A,地点名,ファイル名,日,時,日時,種名,ID_B,雌雄,行動

1,A,A1.avi,2014/01/01,10:00,2014/01/01 10:00,エゾシカ,3,メス,採餌中

2,A,A2.avi,2014/01/01,10:05,2014/01/01 10:05,エゾシカ,3,メス,採餌中

3,A,A3.avi,2014/01/01,11:00,2014/01/01 11:00,エゾタヌキ,3,オス,採餌中

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

ID_A,地点名,ファイル名,日,時,日時,種名,ID_B,雌雄,行動,gid
1,A,A1.avi,2014/01/01,10:00,2014/01/01 10:00,エゾシカ,3,メス,採餌中,1
2,A,A2.avi,2014/01/01,10:05,2014/01/01 10:05,エゾシカ,3,メス,採餌中,1
3,A,A3.avi,2014/01/01,11:00,2014/01/01 11:00,エゾタヌキ,3,オス,採餌中,3

ID_A,地点名,ファイル名,日,時,日時,種名,ID_B,雌雄,行動,gid

1,A,A1.avi,2014/01/01,10:00,2014/01/01 10:00,エゾシカ,3,メス,採餌中,1

2,A,A2.avi,2014/01/01,10:05,2014/01/01 10:05,エゾシカ,3,メス,採餌中,1

3,A,A3.avi,2014/01/01,11:00,2014/01/01 11:00,エゾタヌキ,3,オス,採餌中,3

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

# -*- coding: utf8 -*-
##　標準モジュールのインポート 
import datetime
import csv
import dateutil.parser
### 入力例
## fpath  "C:\\temp\\sample.csv"
## intval  10
def gdata(fpath, intval):
    # csv ファイルの読み取り
    cobj = csv.reader(open(fpath, "rb"))
    # csv ファイルを2次元の配列データに変換
    arry = [ v for v in cobj]
    # 撮影間隔をdatetime形式に変換
    td = datetime.timedelta(0,0,0,0,intval,0,0)
    # 配列データの1行目にgid列を追加   
    ini = arry[0]
    ini.append("gid")
    # 繰り返し処理の初期値(0は列名行のため1からはじめる)
    i = 1
    x = 1
    # レコードの末端まで繰り返し処理
    while i < len(arry):
        # 配列のi行目にiというgidを挿入
        a = arry[i]
        a.append(str(i)) 
        # 次の行(レコード)がなければここで処理終了
        if x + 1 == len(arry):
            break
        # 次の行があれば処理を続行
        x2 = x + 1
        # x行目とx+1行目の6列目のdatetime文字列をdatetime形式にパース
        dto1 = dateutil.parser.parse(arry[x][5])
        dto2 = dateutil.parser.parse(arry[x2][5])
        # 同種同性で前のものと10分以内に撮影された個体であればiというグループ化ID(gid)を付与する
        while dto2 - dto1 < td and arry[x2][6]==arry[x][6] and arry[x2][8]==arry[x][8]:
            x = x + 1
            a = arry[x]
            a.append(str(i))
            # 次の行(レコード)がなければここで処理終了
            if x + 1 == len(arry):
                break
            # 次の行があれば処理を続行
            x2 = x + 1
            # x行目とx+1行目の6列目のdatetime文字列をdatetime形式にパース
            dto1 = dateutil.parser.parse(arry[x][5])
            dto2 = dateutil.parser.parse(arry[x2][5])            
        # 次に作成されるグループのリスト内での順番を定義
        i = x + 1
        x = x + 1
    # 同じフォルダに _g をつけたCSVファイルを新規作成 (同じ名前のファイルがあれば上書きする)
    f = open(fpath[0:-4] + "_g.csv", 'wb')
    # 配列データをCSVファイルに変換
    csvWriter = csv.writer(f)
    csvWriter.writerows(arry)    
    # CSVファイル書き込み終了
    f.close()
    print "finish" 

if __name__ == '__main__':
    print "this is code block"

# -*- coding: utf8 -*-

##　標準モジュールのインポート

import datetime

import csv

import dateutil.parser

### 入力例

## fpath "C:\\temp\\sample.csv"

## intval 10

def gdata(fpath, intval):

# csv ファイルの読み取り

cobj = csv.reader(open(fpath, "rb"))

# csv ファイルを2次元の配列データに変換

arry = [ v for v in cobj]

# 撮影間隔をdatetime形式に変換

td = datetime.timedelta(0,0,0,0,intval,0,0)

# 配列データの1行目にgid列を追加

ini = arry[0]

ini.append("gid")

# 繰り返し処理の初期値(0は列名行のため1からはじめる)

i = 1

x = 1

# レコードの末端まで繰り返し処理

while i < len(arry):

# 配列のi行目にiというgidを挿入

a = arry[i]

a.append(str(i))

# 次の行(レコード)がなければここで処理終了

if x + 1 == len(arry):

break

# 次の行があれば処理を続行

x2 = x + 1

# x行目とx+1行目の6列目のdatetime文字列をdatetime形式にパース

dto1 = dateutil.parser.parse(arry[x][5])

dto2 = dateutil.parser.parse(arry[x2][5])

# 同種同性で前のものと10分以内に撮影された個体であればiというグループ化ID(gid)を付与する

while dto2 - dto1 < td and arry[x2][6]==arry[x][6] and arry[x2][8]==arry[x][8]:

x = x + 1

a = arry[x]

a.append(str(i))

# 次の行(レコード)がなければここで処理終了

if x + 1 == len(arry):

break

# 次の行があれば処理を続行

x2 = x + 1

# x行目とx+1行目の6列目のdatetime文字列をdatetime形式にパース

dto1 = dateutil.parser.parse(arry[x][5])

dto2 = dateutil.parser.parse(arry[x2][5])

# 次に作成されるグループのリスト内での順番を定義

i = x + 1

x = x + 1

# 同じフォルダに _g をつけたCSVファイルを新規作成 (同じ名前のファイルがあれば上書きする)

f = open(fpath[0:-4] + "_g.csv", 'wb')

# 配列データをCSVファイルに変換

csvWriter = csv.writer(f)

csvWriter.writerows(arry)

# CSVファイル書き込み終了

f.close()

print "finish"

if __name__ == '__main__':

print "this is code block"

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

投稿日時: 2014年9月21日投稿者: admin

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)の画像データを撮影間隔でグループ化

投稿日時: 2014年9月15日投稿者: admin

返信

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

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

# -*- coding: utf8 -*-
##　標準モジュールのインポート 
import os
import glob
import datetime
import shutil

### 入力例
## path  "C:\\data\\nakates\\"
## ext  "jpg"
## intval  10
def gpic(path, ext, intval):
    # 対象フォルダ内のファイルリスト作成
    flist = glob.glob(path + '*.' + str(ext))
    # リストの要素の個数をカウント
    fcount = len(flist)
    # リストをソート
    flist.sort()
    # 変数の初期値を設定
    # i はグループ化の名称(グループ内で最も早く撮影されたデータの順序)
    i = 0
    # x はリスト内の順を示す変数
    x = 0
    while i < fcount:
        # ファイルの保存フォルダを作成
        os.mkdir(path + str(i))
        # x番目のファイルを保存フォルダ内にコピー
        shutil.copy2(flist[x],path + str(i) + "\\" + os.path.basename(flist[x]))
        # リストの要素数を超えないように条件設定
        if x + 1 == fcount:
            break
        # x+1番目のファイルをx2と定義
        x2 = x + 1
        # x番目およびx+1番目のファイルのタイムスタンプをdatetime形式で抽出
        dt = datetime.datetime.fromtimestamp(os.stat(flist[x]).st_mtime)
        dt2 = datetime.datetime.fromtimestamp(os.stat(flist[x2]).st_mtime)
        # x番目およびx+1番目のファイルのタイムスタンプの差が10分を超えなければループ
        while dt2 - dt < datetime.timedelta(0,0,0,0,intval,0,0):
            # x の10分後以内にx+1番目のファイルが撮影された場
            # x+1番目のファイルを保存フォルダ内にコピー
            x = x + 1
            shutil.copy2(flist[x],path + str(i) + "\\" + os.path.basename(flist[x]))
            ## x + 1 がリストの要素数と同じ場合ループ終了
            if x + 1 == fcount:
                break
            # x+1番目のファイルをx2と定義
            x2 = x + 1
            # x番目およびx+1番目のファイルのタイムスタンプをdatetime形式で抽出
            dt = datetime.datetime.fromtimestamp(os.stat(flist[x]).st_mtime)
            dt2 = datetime.datetime.fromtimestamp(os.stat(flist[x2]).st_mtime)
        # チェック用メッセージ
        print "grouping segment " + str(i) + " was finished"
        print "group " + str(i) + " is including files from no " + str(i) + " to " + str(x)
        # ディレクトリ内ファイル数をカウント
        pcount = int(x) - int(i) + 1
        # ディレクトリ名をリネーム 　start_end_count ex 15_20_6
        os.rename(path + str(i),path + str(i) + "_" + str(x) + "_" + str(pcount))
        # 次に作成されるグループのリスト内での順番を定義
        i = x + 1
        x = x + 1
    # ディレクトリ名をリネーム(最後のファイル用) 
    os.rename(path + str(i),path + str(i) + "_" + str(x) + "_" + str(pcount))    
    print "finish" 

if __name__ == '__main__':
    print "this is code block"

# -*- coding: utf8 -*-

##　標準モジュールのインポート

import os

import glob

import datetime

import shutil

### 入力例

## path "C:\\data\\nakates\\"

## ext "jpg"

## intval 10

def gpic(path, ext, intval):

# 対象フォルダ内のファイルリスト作成

flist = glob.glob(path + '*.' + str(ext))

# リストの要素の個数をカウント

fcount = len(flist)

# リストをソート

flist.sort()

# 変数の初期値を設定

# i はグループ化の名称(グループ内で最も早く撮影されたデータの順序)

i = 0

# x はリスト内の順を示す変数

x = 0

while i < fcount:

# ファイルの保存フォルダを作成

os.mkdir(path + str(i))

# x番目のファイルを保存フォルダ内にコピー

shutil.copy2(flist[x],path + str(i) + "\\" + os.path.basename(flist[x]))

# リストの要素数を超えないように条件設定

if x + 1 == fcount:

break

# x+1番目のファイルをx2と定義

x2 = x + 1

# x番目およびx+1番目のファイルのタイムスタンプをdatetime形式で抽出

dt = datetime.datetime.fromtimestamp(os.stat(flist[x]).st_mtime)

dt2 = datetime.datetime.fromtimestamp(os.stat(flist[x2]).st_mtime)

# x番目およびx+1番目のファイルのタイムスタンプの差が10分を超えなければループ

while dt2 - dt < datetime.timedelta(0,0,0,0,intval,0,0):

# x の10分後以内にx+1番目のファイルが撮影された場

# x+1番目のファイルを保存フォルダ内にコピー

x = x + 1

shutil.copy2(flist[x],path + str(i) + "\\" + os.path.basename(flist[x]))

## x + 1 がリストの要素数と同じ場合ループ終了

if x + 1 == fcount:

break

# x+1番目のファイルをx2と定義

x2 = x + 1

# x番目およびx+1番目のファイルのタイムスタンプをdatetime形式で抽出

dt = datetime.datetime.fromtimestamp(os.stat(flist[x]).st_mtime)

dt2 = datetime.datetime.fromtimestamp(os.stat(flist[x2]).st_mtime)

# チェック用メッセージ

print "grouping segment " + str(i) + " was finished"

print "group " + str(i) + " is including files from no " + str(i) + " to " + str(x)

# ディレクトリ内ファイル数をカウント

pcount = int(x) - int(i) + 1

# ディレクトリ名をリネーム　start_end_count ex 15_20_6

os.rename(path + str(i),path + str(i) + "_" + str(x) + "_" + str(pcount))

# 次に作成されるグループのリスト内での順番を定義

i = x + 1

x = x + 1

# ディレクトリ名をリネーム(最後のファイル用)

os.rename(path + str(i),path + str(i) + "_" + str(x) + "_" + str(pcount))

print "finish"

if __name__ == '__main__':

print "this is code block"

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

それで

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

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

# モジュールのインポート
from camtrap import gpic as g
# g.gpic("dataディレクトリ","対象ファイルの拡張子",指定する間隔(分単位))
# 入力例
g.gpic("C:\\data\\","jpg",10)

# モジュールのインポート

from camtrap import gpic as g

# g.gpic("dataディレクトリ","対象ファイルの拡張子",指定する間隔(分単位))

# 入力例

g.gpic("C:\\data\\","jpg",10)

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

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

TOPEXをGRASSで算出してみる

投稿日時: 2014年8月12日投稿者: admin

返信

もくもく単純データ処理に嫌気がさしてきたので
気分転換を兼ねて,
以前相談を受けていて放置していた
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だと思ってください…

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

とりあえず

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

という前提で進みます。

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

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

それで

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

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

# モジュールのインポート
from topog import topex as t
# t.topex("gisdataディレクトリ","ロケーション","マップセット",
#　　　　　　"入力DEM名",DEMの解像度(m),計測間隔(m),計測範囲(m),"出力ラスタ名")
# 入力例
t.topex("C:\\temp\\grass","hokk","PERMANENT","dem",10,100,2000,"topex_raster")

# モジュールのインポート

from topog import topex as t

# t.topex("gisdataディレクトリ","ロケーション","マップセット",

#　　　　　　"入力DEM名",DEMの解像度(m),計測間隔(m),計測範囲(m),"出力ラスタ名")

# 入力例

t.topex("C:\\temp\\grass","hokk","PERMANENT","dem",10,100,2000,"topex_raster")

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

## import module ##
import grass.script as grass
import grass.script.setup as gsetup
import math
import os

def topex(gisdb,location,mapset,dem,res,interval,range_max,output):

    # set up grass path #
    gisbase = os.environ['GISBASE']
    # set GRASS gisdatabase, location, mapset #
    gsetup.init(gisbase, gisdb, location, mapset)
    # set region as same as input DEM (extent and resolution)
    grass.run_command("g.region", rast=dem)
    # calculate pixel based interval
    interval_p = int(round(interval / res))
    # calculate pixel based maximum range
    range_max_p = int(round(range_max / res))
    
    ## set formula for calculation ##
    # set initial part of formula
    formula_t1 = output + " = "
    # 8 direction as azimuth angle
    for azi in (0,45,90,135,180,225,270,315):
        # initial cell number
        num = 0
        formula = "max("
        while num &lt; range_max_p:
            num = num + interval_p
            x = int(round(math.sin(math.radians(azi))*num))
            y = int(round(math.cos(math.radians(azi))*num))
            formula = formula + "atan((" + dem + "[" + str(y) + "," + str(x) + "] - " + dem + ") / (" + str(res) + " * " + str(num) + ")),"
        # build formula for calculating max elevation angle at each direction    
        formula_t = formula[0:-1] + ")"
        # build formula for summing up value at all direction 
        formula_t1 = formula_t1 + formula_t + " + "
    # build formula for r.mapcalc
    formula_f = formula_t1[0:-3]    
    print formula_f    
    ## run r.mapcalc command ##
    grass.run_command("r.mapcalc_wrapper",expression = formula_f)
    
    print "finish"
    
if __name__ == '__main__':
    print "this is code block"

## import module ##

import grass.script as grass

import grass.script.setup as gsetup

import math

import os

def topex(gisdb,location,mapset,dem,res,interval,range_max,output):

# set up grass path #

gisbase = os.environ['GISBASE']

# set GRASS gisdatabase, location, mapset #

gsetup.init(gisbase, gisdb, location, mapset)

# set region as same as input DEM (extent and resolution)

grass.run_command("g.region", rast=dem)

# calculate pixel based interval

interval_p = int(round(interval / res))

# calculate pixel based maximum range

range_max_p = int(round(range_max / res))

## set formula for calculation ##

# set initial part of formula

formula_t1 = output + " = "

# 8 direction as azimuth angle

for azi in (0,45,90,135,180,225,270,315):

# initial cell number

num = 0

formula = "max("

while num < range_max_p:

num = num + interval_p

x = int(round(math.sin(math.radians(azi))*num))

y = int(round(math.cos(math.radians(azi))*num))

formula = formula + "atan((" + dem + "[" + str(y) + "," + str(x) + "] - " + dem + ") / (" + str(res) + " * " + str(num) + ")),"

# build formula for calculating max elevation angle at each direction

formula_t = formula[0:-1] + ")"

# build formula for summing up value at all direction

formula_t1 = formula_t1 + formula_t + " + "

# build formula for r.mapcalc

formula_f = formula_t1[0:-3]

print formula_f

## run r.mapcalc command ##

grass.run_command("r.mapcalc_wrapper",expression = formula_f)

print "finish"

if __name__ == '__main__':

print "this is code block"

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

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

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

# モジュールのインポート
from topog import topexf as t
# t.topexf("入力DEM名",DEMの解像度(m),計測間隔(m),計測範囲(m),"出力ラスタ名","出力テキストファイル(フルパス)")
# 入力例
t.topex("dem",10,100,2000,"topex_raster01","C:\\temp\\topex_command.txt")

# モジュールのインポート

from topog import topexf as t

# t.topexf("入力DEM名",DEMの解像度(m),計測間隔(m),計測範囲(m),"出力ラスタ名","出力テキストファイル(フルパス)")

# 入力例

t.topex("dem",10,100,2000,"topex_raster01","C:\\temp\\topex_command.txt")

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

## import math module ##
import math

def topexf(dem,res,interval,range_max,output,out_formula):

    # calculate pixel based interval
    interval_p = int(round(interval / res))
    # calculate pixel based maximum range
    range_max_p = int(round(range_max / res))
    
    ## set formula for calculation ##
    # set initial formula
    formula_t1 = output + " = "
    
    for azi in (0,45,90,135,180,225,270,315):
        # initial cell number
        num = 0
        formula = "max("
        while num &lt; range_max_p:
            num = num + interval_p
            x = int(round(math.sin(math.radians(azi))*num))
            y = int(round(math.cos(math.radians(azi))*num))
            formula = formula + "atan((" + dem + "[" + str(y) + "," + str(x) + "] - " + dem + ") / (" + str(res) + " * " + str(num) + ")),"
        # build formula for calculating max elevation angle at single direction    
        formula_t = formula[0:-1] + ")"
        # build formula for summing up value at all direction 
        formula_t1 = formula_t1 + formula_t + " + "
    # build formula for r.mapcalc
    formula_f = formula_t1[0:-3]    
    print formula_f    

    f = open(out_formula, 'w')
    f.write("r.mapcalc '" + formula_f + "'")
    f.close()
    
    print "finish"
    
if __name__ == '__main__':
    print "this is code block"

## import math module ##

import math

def topexf(dem,res,interval,range_max,output,out_formula):

# calculate pixel based interval

interval_p = int(round(interval / res))

# calculate pixel based maximum range

range_max_p = int(round(range_max / res))

## set formula for calculation ##

# set initial formula

formula_t1 = output + " = "

for azi in (0,45,90,135,180,225,270,315):

# initial cell number

num = 0

formula = "max("

while num < range_max_p:

num = num + interval_p

x = int(round(math.sin(math.radians(azi))*num))

y = int(round(math.cos(math.radians(azi))*num))

formula = formula + "atan((" + dem + "[" + str(y) + "," + str(x) + "] - " + dem + ") / (" + str(res) + " * " + str(num) + ")),"

# build formula for calculating max elevation angle at single direction

formula_t = formula[0:-1] + ")"

# build formula for summing up value at all direction

formula_t1 = formula_t1 + formula_t + " + "

# build formula for r.mapcalc

formula_f = formula_t1[0:-3]

print formula_f

f = open(out_formula, 'w')

f.write("r.mapcalc '" + formula_f + "'")

f.close()

print "finish"

if __name__ == '__main__':

print "this is code block"

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

ではでは。

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

投稿日時: 2014年3月11日投稿者: admin

返信

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

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

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

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

import shapefile
import os
import fnmatch

# データが保存されているパスをセット
dir_p = "C:¥¥temp_ctr¥¥"
# 拡張子を指定(*.shp)
ext = "*.shp"
# *.shpのファイル名を抽出しリストにする
shplist = fnmatch.filter(os.listdir(dir_p), ext)

# データのサブセット処理
try:
    for i in shplist:        
        # データの読み込み
        r = shapefile.Reader(dir_p + str(i))
        # 書き込みインスタンスを定義(ポリライン)
        w = shapefile.Writer(shapeType=shapefile.POLYLINE)
        # 属性フィールドをコピー
        w.fields = list(r.fields)
        # 空のリストをつくる 
        sel_list = [] 
        for rec in enumerate(r.records()):
            # 100の倍数の標高値(8列目の属性フィールドに格納されている)を持つ等高線を抽出し書き込む
            if rec[1][7] % 100 == 0:
                sel_list.append(rec) 
        # ジオメトリと属性値を書き込む
        for rec in sel_list:
            w._shapes.append(r.shape(rec[0]))
            w.records.append(rec[1])
        # shapefileとして保存する(ファイル名の先頭にctr100mを付与)
        w.save(dir_p + "ctr100m" + str(i))
except:
    print "error!"
print "finish all"

import shapefile

import os

import fnmatch

# データが保存されているパスをセット

dir_p = "C:¥¥temp_ctr¥¥"

# 拡張子を指定(*.shp)

ext = "*.shp"

# *.shpのファイル名を抽出しリストにする

shplist = fnmatch.filter(os.listdir(dir_p), ext)

# データのサブセット処理

try:

for i in shplist:

# データの読み込み

r = shapefile.Reader(dir_p + str(i))

# 書き込みインスタンスを定義(ポリライン)

w = shapefile.Writer(shapeType=shapefile.POLYLINE)

# 属性フィールドをコピー

w.fields = list(r.fields)

# 空のリストをつくる

sel_list = []

for rec in enumerate(r.records()):

# 100の倍数の標高値(8列目の属性フィールドに格納されている)を持つ等高線を抽出し書き込む

if rec[1][7] % 100 == 0:

sel_list.append(rec)

# ジオメトリと属性値を書き込む

for rec in sel_list:

w._shapes.append(r.shape(rec[0]))

w.records.append(rec[1])

# shapefileとして保存する(ファイル名の先頭にctr100mを付与)

w.save(dir_p + "ctr100m" + str(i))

except:

print "error!"

print "finish all"

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

投稿日時: 2013年12月2日投稿者: admin

返信

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

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

>>FOSS4G 勉強会 in 函館

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

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

これは 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も変更してください。

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

以下より

変数名変数値

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

GISBASE C:¥OSGeo4W¥apps¥grass¥grass-6.4.3
GISRC C:¥Users¥username¥AppData¥Roaming¥GRASS6¥grass6
LD_LIBRARY_PATH C:¥OSGeo4W¥apps¥grass¥grass-6.4.3¥lib

PATH
C:¥OSGeo4W¥apps¥grass¥grass-6.4.3¥etc
C:¥OSGeo4W¥apps¥grass¥grass-6.4.3¥etc¥python
C:¥OSGeo4W¥apps¥grass¥grass-6.4.3¥lib
C:¥OSGeo4W¥apps¥grass¥grass-6.4.3¥bin
C:¥OSGeo4W¥apps¥grass¥grass-6.4.3¥extralib
C:¥OSGeo4W¥apps¥grass¥grass-6.4.3¥msys¥bin

PYTHONLIB C:¥Python27
PYTHONPATH C:¥OSGeo4W¥apps¥grass¥grass-6.4.3¥etc¥python
GRASS_SH C:¥OSGeo4W¥apps¥msys¥bin¥sh.exe

GISBASE C:¥OSGeo4W¥apps¥grass¥grass-6.4.3

GISRC C:¥Users¥username¥AppData¥Roaming¥GRASS6¥grass6

LD_LIBRARY_PATH C:¥OSGeo4W¥apps¥grass¥grass-6.4.3¥lib

PATH

C:¥OSGeo4W¥apps¥grass¥grass-6.4.3¥etc

C:¥OSGeo4W¥apps¥grass¥grass-6.4.3¥etc¥python

C:¥OSGeo4W¥apps¥grass¥grass-6.4.3¥lib

C:¥OSGeo4W¥apps¥grass¥grass-6.4.3¥bin

C:¥OSGeo4W¥apps¥grass¥grass-6.4.3¥extralib

C:¥OSGeo4W¥apps¥grass¥grass-6.4.3¥msys¥bin

PYTHONLIB C:¥Python27

PYTHONPATH C:¥OSGeo4W¥apps¥grass¥grass-6.4.3¥etc¥python

GRASS_SH C:¥OSGeo4W¥apps¥msys¥bin¥sh.exe

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

>>スクリプトを書くよ

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

## モジュールのインポート

import grass.script as grass
import grass.script.setup as gsetup
import os

## GRASSのデータベースのセットアップ
# GRASSのプログラムが格納されているディレクトリを指定

gisbase = os.environ['GISBASE']

# GRASSのgisdatabase, location, mapsetを指定
# *GUIで使用する際に最初に指定する項目と一緒です

gisdb = "C:/foss4ghakodate/grass"
location = "hakodate"
mapset = "srtm"
gsetup.init(gisbase, gisdb, location, mapset)

## モジュールのインポート

import grass.script as grass

import grass.script.setup as gsetup

import os

## GRASSのデータベースのセットアップ

# GRASSのプログラムが格納されているディレクトリを指定

gisbase = os.environ['GISBASE']

# GRASSのgisdatabase, location, mapsetを指定

# *GUIで使用する際に最初に指定する項目と一緒です

gisdb = "C:/foss4ghakodate/grass"

location = "hakodate"

mapset = "srtm"

gsetup.init(gisbase, gisdb, location, mapset)

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

## GRASSのコマンドを用いた処理

# 変数を定義

srtm = "hakodate_j2ku54N"
stat = "variance"
expdir = "C:/temp/"

# 以下の処理を実行
# r.neighborsで近傍解析を行い
# r.out.gdal でgeotif形式でエクスポート

for diam in(3,5,7,11,13,15):
    outdata =  "focal_dem_" + str(diam) + "x" + str(diam) + "_" + stat[0:3]
    expdata = expdir + outdata + ".tif"
    grass.run_command("r.neighbors", input = srtm, output = outdata, method = stat, size = diam)
    grass.run_command("r.out.gdal", input = outdata, output = expdata)
    print outdata

print "finish"

## GRASSのコマンドを用いた処理

# 変数を定義

srtm = "hakodate_j2ku54N"

stat = "variance"

expdir = "C:/temp/"

# 以下の処理を実行

# r.neighborsで近傍解析を行い

# r.out.gdal でgeotif形式でエクスポート

for diam in(3,5,7,11,13,15):

outdata = "focal_dem_" + str(diam) + "x" + str(diam) + "_" + stat[0:3]

expdata = expdir + outdata + ".tif"

grass.run_command("r.neighbors", input = srtm, output = outdata, method = stat, size = diam)

grass.run_command("r.out.gdal", input = outdata, output = expdata)

print outdata

print "finish"

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

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

投稿日時: 2012年12月13日投稿者: admin

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

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

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

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

以降

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

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

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

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

webtisMap = new webtis.Layer.BaseMap("webtismap");
map.addLayer(webtisMap);

1 2	webtisMap = new webtis.Layer.BaseMap("webtismap"); map.addLayer(webtisMap);

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

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

webtisMap = new webtis.Layer.BaseMap("webtismap");
var dataSet = webtisMap.getOrthoDataSet();
webtisMap.setDataSet(dataSet);
map.addLayer(webtisMap);

webtisMap = new webtis.Layer.BaseMap("webtismap");

var dataSet = webtisMap.getOrthoDataSet();

webtisMap.setDataSet(dataSet);

map.addLayer(webtisMap);

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

そして

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

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

        webtisMap = new webtis.Layer.BaseMap("GSI CyberJapanMap");
        map.addLayer(webtisMap);

1 2	webtisMap = new webtis.Layer.BaseMap("GSI CyberJapanMap"); map.addLayer(webtisMap);

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

        webtisMap = new webtis.Layer.BaseMap("GSI CyberJapanMap");
	var dataSet = webtisMap.getOrthoDataSet();
	webtisMap.setDataSet(dataSet);
        map.addLayer(webtisMap);

webtisMap = new webtis.Layer.BaseMap("GSI CyberJapanMap");

var dataSet = webtisMap.getOrthoDataSet();

webtisMap.setDataSet(dataSet);

map.addLayer(webtisMap);

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の中身は以下の通り

&lt;html xmlns="http://www.w3.org/1999/xhtml"&gt;
  &lt;head&gt;
    &lt;title&gt;OpenLayers GSI CyberJapanMap Layer&lt;/title&gt;
    &lt;link rel="stylesheet" href="qgis.css" type="text/css"&gt;
    &lt;script src="http://portal.cyberjapan.jp/sys/OpenLayers-2.11/OpenLayers.js"&gt;&lt;/script&gt;
    &lt;script src="http://portal.cyberjapan.jp/sys/v4/webtis/webtis_v4.js"&gt;&lt;/script&gt;
    &lt;script src="OlOverviewMarker.js"&gt;&lt;/script&gt;
    &lt;script type="text/javascript"&gt;
      var map;
      var loadEnd;
      var oloMarker; // OpenLayer Overview Marker
      function init() {
        map = new OpenLayers.Map('map', {
          theme: null,
          controls: [],
          projection: new OpenLayers.Projection("EPSG:3857"),
        });

        loadEnd = false;
        function layerLoadStart(event) {
          loadEnd = false;
        }
        function layerLoadEnd(event) {
          loadEnd = true;
        }

        webtisMap = new webtis.Layer.BaseMap("GSI CyberJapanMap");
	var dataSet = webtisMap.getOrthoDataSet();
	webtisMap.setDataSet(dataSet);
        map.addLayer(webtisMap);

        map.addControl(new OpenLayers.Control.Attribution());
        oloMarker = new OlOverviewMarker(map, getPathUpper(document.URL) + '/x.png');
      }
    &lt;/script&gt;
  &lt;/head&gt;
  &lt;body onload="init()"&gt;
    &lt;div id="map"&gt;&lt;/div&gt;
  &lt;/body&gt;
&lt;/html&gt;

<head>

<title>OpenLayers GSI CyberJapanMap Layer</title>

var map;

var loadEnd;

var oloMarker; // OpenLayer Overview Marker

function init() {

map = new OpenLayers.Map('map', {

theme: null,

controls: [],

projection: new OpenLayers.Projection("EPSG:3857"),

});

loadEnd = false;

function layerLoadStart(event) {

loadEnd = false;

}

function layerLoadEnd(event) {

loadEnd = true;

}

webtisMap = new webtis.Layer.BaseMap("GSI CyberJapanMap");

var dataSet = webtisMap.getOrthoDataSet();

webtisMap.setDataSet(dataSet);

map.addLayer(webtisMap);

map.addControl(new OpenLayers.Control.Attribution());

oloMarker = new OlOverviewMarker(map, getPathUpper(document.URL) + '/x.png');

}

</script>

</head>

</body>

</html>

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

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

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

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

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

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

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

&lt;html xmlns="http://www.w3.org/1999/xhtml"&gt;
  &lt;head&gt;
    &lt;title&gt;OpenLayers GSI CyberJapanMap Layer&lt;/title&gt;
    &lt;link rel="stylesheet" href="qgis.css" type="text/css"&gt;
    &lt;script src="http://portal.cyberjapan.jp/sys/OpenLayers-2.11/OpenLayers.js"&gt;&lt;/script&gt;
    &lt;script src="http://portal.cyberjapan.jp/sys/v4/webtis/webtis_v4.js"&gt;&lt;/script&gt;
    &lt;script src="OlOverviewMarker.js"&gt;&lt;/script&gt;
    &lt;script type="text/javascript"&gt;
      var map;
      var loadEnd;
      var oloMarker; // OpenLayer Overview Marker
      function init() {
        map = new OpenLayers.Map('map', {
          theme: null,
          controls: [],
          projection: new OpenLayers.Projection("EPSG:3857"),
        });

        loadEnd = false;
        function layerLoadStart(event) {
          loadEnd = false;
        }
        function layerLoadEnd(event) {
          loadEnd = true;
        }

	var dataSet = {
		        5:{dataId:"JAIS"},
		        6:{dataId:"JAIS"},
		        7:{dataId:"JAIS"},
		        8:{dataId:"JAIS"},
		        9:{dataId:"BAFD1000K"},
		        10:{dataId:"BAFD1000K"},
		        11:{dataId:"BAFD1000K"},
		        12:{dataId:"BAFD200K"},
		        13:{dataId:"BAFD200K"},
		        14:{dataId:"BAFD200K"},
		        15:{dataId:"NLII1"},
		        16:{dataId:"NLII1"},
		        17:{dataId:"NLII1"} 
        	};

        webtisMap = new webtis.Layer.BaseMap("GSI CyberJapanMap", {
	dataSet:dataSet
	});

        map.addLayer(webtisMap);

        map.addControl(new OpenLayers.Control.Attribution());
        oloMarker = new OlOverviewMarker(map, getPathUpper(document.URL) + '/x.png');
      }
    &lt;/script&gt;
  &lt;/head&gt;
  &lt;body onload="init()"&gt;
    &lt;div id="map"&gt;&lt;/div&gt;
  &lt;/body&gt;
&lt;/html&gt;

<head>

<title>OpenLayers GSI CyberJapanMap Layer</title>

var map;

var loadEnd;

var oloMarker; // OpenLayer Overview Marker

function init() {

map = new OpenLayers.Map('map', {

theme: null,

controls: [],

projection: new OpenLayers.Projection("EPSG:3857"),

});

loadEnd = false;

function layerLoadStart(event) {

loadEnd = false;

}

function layerLoadEnd(event) {

loadEnd = true;

}

var dataSet = {

5:{dataId:"JAIS"},

6:{dataId:"JAIS"},

7:{dataId:"JAIS"},

8:{dataId:"JAIS"},

9:{dataId:"BAFD1000K"},

10:{dataId:"BAFD1000K"},

11:{dataId:"BAFD1000K"},

12:{dataId:"BAFD200K"},

13:{dataId:"BAFD200K"},

14:{dataId:"BAFD200K"},

15:{dataId:"NLII1"},

16:{dataId:"NLII1"},

17:{dataId:"NLII1"}

};

webtisMap = new webtis.Layer.BaseMap("GSI CyberJapanMap", {

dataSet:dataSet

});

map.addLayer(webtisMap);

map.addControl(new OpenLayers.Control.Attribution());

oloMarker = new OlOverviewMarker(map, getPathUpper(document.URL) + '/x.png');

}

</script>

</head>

</body>

</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行を挿入

    #　電子国土基本図（オルソ画像）表示用
    self.olLayerTypeRegistry.add( OlLayerType(self, 'GSI CyberJapanOrtho', 'cjp_icon.png', 'cjp_ortho.html', True) )

    #　国土画像情報第1期(1974年～1978年)表示用
    self.olLayerTypeRegistry.add( OlLayerType(self, 'GSI CyberJapanOrtho70s', 'cjp_icon.png', 'cjp_ortho70s.html', True) )

#　電子国土基本図（オルソ画像）表示用