GRASSを使ってみる

post2008,03,03 00:56 GMT

QGIS for Windows

三重大学生物資源学部の地理情報システム学２にQGIS for Windowsの紹介がありました。

試しにVMWare上のWindows XPにQuantum GIS 0.9 Downloadsからダウンロードした0.9.2rc1をインストールしてみたところ、非常に楽にGRASS環境を導入できました。。

(08/03/05 追記 まだ充分、環境が整っていないようです。しばらく、最新版を追いかければ使える環境になるでしょうか。具体的には、Grass toolsのモジュールで動かないものがいくつもあることや、コマンドラインできっちり動作させられないコマンドなどがあります。)

linux上のQGISでも同様な作業ができましたので、そこからGRASSを操作する方法で、石垣島のデータセットを整える例をお見せします。

qgis自体のインストールをc:\qgisとしました。

qgisを動かしているベースシステムであるmsysでは、コンソールからアクセスできるディレクトリがc:\qgis\msys\以下に限られ ますので、grassにインポートするデータなどはc:\qgis\msys\home\”user name”以下に一時的に保存します。

その後、データセットを保存していくためのディレクトリを決めます。

c:\にgis-dataというディレクトリを作成します。UNIX系のソフトを使用する際には、ディレクトリ名にスペースが入ったり、アスキー文字以外が使われていると思わぬエラーを招きますので、ドキュメントフォルダ以外に用意しておいたほうがいいでしょう。

キャプチャ画面はlinux上でホームディレクトリ以下にgis-dataを作っていますので、読み替えて見てください。

qgisを起動します。

読み込まれているプラグインが不十分なので、ロードします。「プラグイン-Plugin Manager…」を選択します。Plugin Managerが起動するので「すべてを選択」ボタンをクリックします。「OK」で、grass pluginがロードされます。

ウィンドウの右端に並んでいるのが、grass pluginのアイコンです。

この時点では、マップセットの選択がされていません。

新しいmapsetを作成します。

New mapsetアイコンをクリックすると、ウィザードが開始されます。

databaseとかかれている項目がGisdbaseになります。ここではc:\gis-dataと入力します。「Next」

Locationの選択になりますが、まだLocationは作られていませんので、「Create New Location」に「Ishigaki-utm」と入力します。

「Next」

Coordinate systemを指定します。「投影法」ボタンを押すと選択肢が選べますので、Projected Coordinate Systemsのツリーを開きます。

ツリーを下にたぐっていくと、Universal Transverse Mercator(UTM)のツリーが出てきます。それを開きます。

かなり、下の方にスクロールするとWGS 84/UTM zone 52Nの選択肢が現れますので選択します。

下部にSRIDやら説明やらが表示されます。

「Next」

Default Regionの指定ダイアログが出てきます。国ごとに指定が可能なので「afghanistan」と表示されたボックスをクリックします。

「Japan」を選択します。「設定」をクリックすると世界地図上に範囲が表示されます。無駄に広い範囲になるのでここで細かく調整してもいいかと思います。

ここでは、そのままで「Next」

mapsetの名前を決めますが、grassでは普通user nameを使います。

「Next」

ここまで指定した形式でディレクトリツリーが作成されます。

「Finish」

作成されました。

これでmapsetのなかに入れましたので、既存のデータをどんどんインポートしていきます。

(08/03/05 追記 上でも記載しましたが、Windows版ではまだ不具合がたくさんあるようです。問題点は追々クリアされるものと思いますが、現時点では以下の作業はlinuxで行っています。)

数値地図(空間データ基盤)の閲覧（試験公開）から数値地図25000をダウンロードします。「閲覧する市区町村の選択」から「沖縄県」をクリックします。

目的の市町村のデータをダウンロードしますが、ここでは石垣市の数値地図25000を選択します。

テンの趣味からsu25000v10.lzhをダウンロードし解凍します。有益なツールに深く感謝します。

su25000.exeを起動し、追加ボタンで47207.lzhを追加します。ArcInfo shapeのラジオボタンをチェックした後実行をするとlzhファイルと同じディレクトリにshpファイルがいくつか作成されています。

47207CM.shp 地名

47207GD.shp 市郡区名

47207GK.shp 確定境界

47207KJ.shp 基準点

47207KK.shp 河川

47207KS.shp

47207SK.shp 水涯線

すべてのデータが変換されているわけではないようです。

「Open GRASS tools」をクリックします。

とてもたくさんのGrassコマンドがモジュールになっています。コマンドラインに慣れていない間はここからコマンドを実行するようにしたらいいですね。

「GRASS shell」を選択すると、linuxではGRASS toolsの中にコンソールのタブができてそこでshellが実行されます。

Comments (0)

post2007,10,12 09:27 IST

DEM Export partII

これまで、GRASSの外にDEMを持ち出す必要がなかったため、DEM Exportを検証できていませんでした。

ところが、location utmの中で5mの解像度という設定で作ったDEMをr.projでlatlonに移動しようと思ったら、移動先の設定解像度に変更されてしまいました。 r.proj不便だぞ。出力側の解像度をキープするオプションがありそうなものだけど。(たぶん、あるんだろうけど)

試行錯誤。

g.region rast=dem5m.ishigaki -p

projection: 1 (UTM)

zone: 52

datum: wgs84

ellipsoid: wgs84

north: 2737340

south: 2698590

west: -906

east: 30179

nsres: 5

ewres: 5

rows: 7750

cols: 6217

cells: 48181750

これを移動したい。

移動先のg.region設定は

projection: 3 (Latitude-Longitude)

zone: 0

datum: ** unknown (default: WGS84) **

ellipsoid: grs80

north: 24:40:20.081262N

south: 24:19:40.080719N

west: 124:04E

east: 124:21:35.000882E

nsres: 0:00:00.081264

ewres: 0:00:00.088887

rows: 15259

cols: 11869

cells: 181109071

この二つはまったく同じ範囲を示しているが、rows、colsをみると解像度がまったく違うことがわかる。

解像度の設定を変更してみる。

g.region res=0.00004712 -p

projection: 3 (Latitude-Longitude)

zone: 0

datum: ** unknown (default: WGS84) **

ellipsoid: grs80

north: 24:40:20.081262N

south: 24:19:40.080719N

west: 124:04E

east: 124:21:35.000882E

nsres: 0:00:00.169631

ewres: 0:00:00.169642

rows: 7310

cols: 6219

cells: 45460890

この様にすると、colsでは近づいたがrowsの方がだいぶ異なる。ということは一画素の表す範囲が違うということになる。これではどうしても細かいところまで詰められない。せっかく5x5mで作ったのに。

試しにr.projをしてみると

r.proj input=dem5m.ishigaki location=wgs mapset=HOGE output=testdem5m.ishigaki

Input Projection Parameters: +proj=utm +zone=52 +a=6378137 +rf=298.257223563 +no_defs +towgs84=0.000,0.000,0.000

Input Unit Factor: 1

Output Projection Parameters: +proj=longlat +a=6378137 +rf=298.257222101 +no_defs

Output Unit Factor: 1

Input:

Cols: 6217 (6217)

Rows: 7750 (7750)

North: 2737340.000000 (2737340.000000)

South: 2698590.000000 (2698590.000000)

West: -906.000000 (-906.000000)

East: 30179.000000 (30179.000000)

EW-res: 5.000000

NS-res: 5.000000

Output:

Cols: 6219 (6219)

Rows: 7310 (7310)

North: 24.672245 (24.672245)

South: 24.327800 (24.327800)

West: 124.066667 (124.066667)

East: 124.359722 (124.359722)

EW-res: 0.000047

NS-res: 0.000047

手作業では精度に疑問が出てしまいます。

そこで、Exporting GRASS Raster Data to ArcGISを参考に

r.out.arc input=dem5m.ishigaki output=dem5m_ishigaki.grd

gdalwarp -co TFW=YES -co INTERLEAVE=PIXEL -s_srs "+proj=utm +zone=52 +a=6378137 +rf=298.257223563 +no_defs +towgs84=0.000,0.000,0.000" -t_srs "+proj=longlat +a=6378137 +rf=298.257222101 +no_defs" -order 1 dem5m_ishigaki.grd dem5m_ishigaki.tif

おおっ、できた。

r.in.gdal input=dem5m_ishigaki.tif output=dem5m.ishigaki

g.region -p rast=dem5m.ishigaki

projection: 3 (Latitude-Longitude)

zone: 0

datum: ** unknown (default: WGS84) **

ellipsoid: grs80

north: 24:40:43.741476N

south: 24:19:13.253866N

west: 124:03:11.921039E

east: 124:22:20.220562E

nsres: 0:00:00.168669

ewres: 0:00:00.168669

rows: 7651

cols: 6808

cells: 52088008

utmからlatlonに投影変換する際に形状が崩れますので、rows colsが違っていますが、こっちの方が信頼できそうですね。

GDAL Utilities

07/10/26 追記

UTMのロケーションで作った地図を、他の環境の人に読める形式で提出したい。

こんな感じでいけますでしょうか。

5mの分解能で渡したい。

g.region -p res=5

projection: 1 (UTM)

zone: 52

datum: wgs84

ellipsoid: wgs84

north: 2730725

south: 2699155

west: -585

east: 28965

nsres: 5

ewres: 5

rows: 6314

cols: 5910

cells: 37315740

d.zoomで範囲を指定。

r.out.arc input=a_map output=a_map.grd

gdalinfo a_map.grd

Driver: AAIGrid/Arc/Info ASCII Grid

Size is 5910, 6314

Coordinate System is `’

Origin = (-585.000000000000000,2730725.000000000000000)

Pixel Size = (5.000000000000000,-5.000000000000000)

Corner Coordinates:

Upper Left ( -585.000, 2730725.000)

Lower Left ( -585.000, 2699155.000)

Upper Right ( 28965.000, 2730725.000)

Lower Right ( 28965.000, 2699155.000)

Center ( 14190.000, 2714940.000)

Band 1 Block=5910x1 Type=Float32, ColorInterp=Undefined

NoData Value=-9999

投影法のデータが含まれていないので、gdalwarpにエクスポート側の投影法を指定した上で、別の投影法に変換する。ここでは、JGD2000のEPSGコード4612を指定している。

gdalwarp -co TFW=YES -co INTERLEAVE=PIXEL

-order 1 -s_srs “`g.proj -wef`” -t_srs EPSG:4612 a_map.grd a_map_jgd.tif

gdalinfo a_map_jgd.tif

Driver: GTiff/GeoTIFF

Size is 6389, 6220

Coordinate System is:

GEOGCS[”JGD2000″,

DATUM[”Japanese_Geodetic_Datum_2000″,

SPHEROID[”GRS 1980″,6378137,298.2572221010042,

AUTHORITY[”EPSG”,”7019″]],

AUTHORITY[”EPSG”,”6612″]],

PRIMEM[”Greenwich”,0],

UNIT[”degree”,0.0174532925199433],

AUTHORITY[”EPSG”,”4612″]]

Origin = (124.058819931968216,24.618900411744811)

Pixel Size = (0.000047165204866,-0.000047165204866)

Metadata:

AREA_OR_POINT=Area

Corner Coordinates:

Upper Left ( 124.0588199, 24.6189004) (124d 3′31.75″E, 24d37′8.04″N)

Lower Left ( 124.0588199, 24.3255328) (124d 3′31.75″E, 24d19′31.92″N)

Upper Right ( 124.3601584, 24.6189004) (124d21′36.57″E, 24d37′8.04″N)

Lower Right ( 124.3601584, 24.3255328) (124d21′36.57″E, 24d19′31.92″N)

Center ( 124.2094892, 24.4722166) (124d12′34.16″E, 24d28′19.98″N)

Band 1 Block=6389x1 Type=Float32, ColorInterp=Gray

これでいいんだけど、r.out.arcでエクスポートしたときにNULLが-9999に変換されちゃうんだよな。他のソフトでは問題にならないんだろうか。

07/10/31 追記

gdalwarpなんかを覚えたら、latlonのロケーションが必要ない感じですね。

latlonで提供されているベースマップも変換して入れちゃえばいいんだし、解析などはUTMでやっちゃおうと思っているので。

gdalinfoだとか、ogrinfoだとか早く覚えておけばよかった。

Comments (0)

post2007,10,01 09:08 IST

New Location

最初の一歩として、新しいlocationを作ってベースマップを整備することをしてみましょう。

07/10/31 修正

途中経過がふっとんでしまっていたので、全面的に書き換えます。前回はWindows版を基にしていましたが、今回はLinux版でいきます。バージョンが違うのでユーザーインターフェイスが異なりますが、自分用のメモということで。

また、latlonのロケーションを作らずutmオンリーでいけるか試してみます。

PERMANENTとuserのmapsetは分けるべきなのでしょうが、あまり理解できていないため全部PERMANENTのmapsetに入れていきます。

ホームディレクトリでgis用のディレクトリを作成します。

mkdir gis

cd gis

grassを起動します。

Define new location with… からEPSGコードを使った作成を選びます。

location nameを指定して、epsgファイルの位置を確認します。

epsgコード表から作りたいロケーションのepsgコードを見つけ出し、ロケーションを作成します。

これで、grassの中に入れます。

okinawaさんがRjpWikiで公開されている、ESRIジャパン株式会社の全国市区町村界データを加工したデータを使います。(著作権表記をご覧ください。)

http://www.okada.jp.org/RWiki/index.php?ShapeFile%A5%E9%A5%A4%A5%D6%A5%E9%A5%EA

の中にある、okinawhk.zipをダウンロードして保存します。

unzipで解凍します。

okinawhk.dbf okinawhk.shp okinawhk.shx

のファイルが現れます。

これは、JGD2000のファイルなのでそのままではインポートできません。utm環境にインポートできるように変換します。

ogr2ogr -s_srs EPSG:4612 -t_srs EPSG:32652 okinawhk_utm.shp okinawhk.shp

s_srsにはこのshpファイルの位置情報がどういう投影法を基にしているのかということを指定します。普通はこの辺の情報を持っているはずな んですが、このファイルには含まれていないので指定する必要がありました。t_srsはターゲットの投影パラメーターです。この場合ロケーションを作成し たときに指定したEPSGコードでいいです。

ここで注意しないといけないのは、このogr2ogrのパラメーターの指定順序は出力ファイルの方が先になるということです。変なの。

うまく変換できたので、インポートしましょう。

v.in.ogr dsn=okinawhk_utm.shp output=okinawa_isl

で入ってくれるはずです。

d.mon x0

g.region vect=okinawa_isl -p

d.rast okinawa_isl

で、確認してみましょう。

地形図もインポートしてみます。

create DEM partIIIの手順通り、Geotiffを作成します。

import-tif.shを少し変更します。上で作成したGeotiffはlatlonなので、utmに変換しながらインポートします。

import-tif.sh

#!/bin/sh

for file in *.tif

do

gdalwarp -s_srs EPSG:4612 -t_srs EPSG:32652 -order 1 $file ${file%.tif}_utm.tif

r.in.gdal -o input=${file%.tif}_utm.tif output=${file%.tif}.m

done

数値地図(空間データ基盤)の閲覧（試験公開）のデータを手元で閲覧できるようにしましょう。

閲覧する市区町村の選択から目的とする地域を含む県を指定します。沖縄県。

市町村のリストが現れますので、必要な市町村を選択します。試しに那覇市。

ファイルサイズのカラムをクリックすると47201.lzhがダウンロードされます。

lha x 47201.lzh

で解凍されます。

できたディレクトリに移動。

変換ソフトを検索します。世の中には素晴らしい先輩方が存在します。Sdf2ogrをビルドしてインストールすると、shapeファイルに変換できます。

データの展開されたディレクトリで

sdf2ogr -s 2457

と実行するとdouro、eki、gaikuなど要素別にディレクトリが掘られ、その中にshape形式のファイルが格納されています。

投影法を確認してみましょう。

ogrinfo -so doro_arc.shp doro_arc

INFO: Open of `doro_arc.shp’

using driver `ESRI Shapefile’ successful.

Layer name: doro_arc

Geometry: Line String

Feature Count: 29094

Extent: (14004.000000, 18436.300000) - (64377.600000, 78673.900000)

Layer SRS WKT:

PROJCS[”JGD2000 / Japan Plane Rectangular CS XV”,

GEOGCS[”JGD2000″,

DATUM[”Mean_Sea_Level”,

SPHEROID[”GRS_1980″,6378137,298.257222101]],

PRIMEM[”Greenwich”,0],

UNIT[”Degree”,0.017453292519943295]],

PROJECTION[”Transverse_Mercator”],

PARAMETER[”latitude_of_origin”,26],

PARAMETER[”central_meridian”,127.5],

PARAMETER[”scale_factor”,0.9999],

PARAMETER[”false_easting”,0],

PARAMETER[”false_northing”,0],

UNIT[”Meter”,1]]

code: String (80.0)

number: String (80.0)

count: Real (11.0)

opt1: String (80.0)

opt2: String (80.0)

opt3: String (80.0)

opt4: String (80.0)

opt5: String (80.0)

やはり、JGD2000です。

じゃあ、変換。

ogr2ogrは変換先のファイルを先に指定する変態的なソフトだということを忘れずに。

ogr2ogr -s_srs EPSG:4612 -t_srs EPSG:32652 doro_arc_utm.shp doro_arc.shp

Comments (0)

post 05:17 IST

tutorial example

大阪市立大学の升本先生、ベンカテッシュ ラガワン先生によるチュートリアルはGRASSを理解するためのものとして有益です。が、GRASS5を対象としているため、現状に合わない部分が出てきています。

スタートでつまづかないように、GRASSの起動の部分だけお見せします。

cygwinGRASSの手順で、GRASS6.1cvsをインストールしてあるものとします。

user nameはHOGEでc:\cygwin\home\HOGEができています。

エクスプローラーでc:\cygwin\home\HOGEにgisというフォルダを作ります。その中に、leics data setをダウンロードし、保存します。

cygwinを起動します。コンソールの中でstartxとします。もう一つターミナル画面が現れます。こちらを使っていきますので、先のコンソールウィンドウは最小化しておきます。

cd gis でディレクトリを移動しておきます。ここにはleics.tar.gzが保存されているはずですので、tar xzvf leics.tar.gzとして解凍します。

grass61と入力し、リターンを押します。コンソールにウェルカム画面が出ます。Enterを入力します。

location、mapsetの選択ダイアログが出ます。

/ 4″ r.mapcalc “b2 = (”mosaic” % 4) / 2″

r.mapcalc “b1 = (”mosaic” % 2) / 1″

b4が等高線情報の入ったレイヤーなのでそれを使いますが、そのままではベクターに出来ません。すべてのラインが1ドットで表現されるように細線化します。

r.thin input=b4 output=contour

それをベクターに変換。

r.to.vect input=contour output=contour

元になった地形図画像と比べてみます。

文字に隠れていた情報も残っているのがわかります。すごいなぁ。

データベースにどのように格納されているのかを確認してみます。

v.info -c contour

INTEGER|cat

INTEGER|value

CHARACTER|label

labelのカラムに標高値を指定していきますが、GRASS6になってv.digitが大幅にリニューアルされたときに便利な機能が取り除かれています。GRASS5のv.digitを使いましょう。

“Some of the useful functions which were available in the first version of the v.digit module has never been implemented in the new versions, e.g. generation of the attribute of isolines (contour height).”

ベクターと参照用の地形図ラスターをエクスポート

GRASS6のベクターはGRASS5と互換性がないため一旦shpファイルにします。このshpファイルをそのまま使うと海岸線のない等高線図になります。海岸線をマージしましょう。

v.out.ogr input=contour dsn=contour.shp format=ESRI_Shapefile

ひぇ〜、4時間ほどかけて80MBのshpファイルが出来上がりました。(ここではかなりの広範囲を扱っています)(07/09/26 追記 ここまでの時点で、一度v.build.polylinesを実行しておくとこんなに時間がかからないかもしれません。)

標高を指定する際、参照できるラスター画像があるとものすごく楽なので、地形図画像もエクスポートします。

g.region rast=mosaic

r.out.tiff input=mosaic output=mosaic.tif -t compression=deflate

• PostGISで

等高線ベクターを海岸線ベクターとマージ

shpファイルをGRASS5に読み込む前に、ここで海岸線を差し込みます。

PostGISを使います。

shp2pgsql -s 4612 -W CP932 -d contour.shp contour>contour.sql

psql -f contour.sql map

qgisでどのように見えるかやってみましょう。

ラスターの色がビビッドでおかしいですが、ちゃんと地図に乗っています。

数値地図25000の水域界.shpが県内の海岸線を含んでいるので、そこから抜き出します。

syuruiのカラムに「水涯線または湖岸線」、「海岸線」、「河口」の3種類が含まれていますが、そのうちの「水涯線または湖岸線」にあたるラインを削除します。

shp2pgsql -s 4612 -W CP932 -d 水域界.shp suiiki > suiiki.sql

psql -f suiiki.sql map

これで、PostGISにsuiikiテーブルが入りました。

pgadmin3などを使って

DELETE FROM suiiki WHERE syurui='水涯線または湖岸線';

海に面したラインが残ります。

ここで地域を指定して、離島の海岸線などを除去してもいいのですが、SQLで指定する方法がわからないため放置しています。後で標高指定をしていないラインをすべて除いてしまいますので、害はありません。GRASS5のv.digitでも取り除けます。

必要のないカラムを削除します。

ALTER TABLE suiiki DROP COLUMN syurui;

ALTER TABLE suiiki DROP COLUMN id;

contourのテーブルでは構造がどうなっているのか確認します。

labelの欄が空ですね。どのラインもまだ標高を指定されていないので当然ですね。そこを使いましょう。

ALTER TABLE suiiki ADD COLUMN label varchar(80);

同じようなカラムを作ります。

これで、gid、label、the_geomのカラムになります。

contourについても同様にします。

ALTER TABLE contour DROP COLUMN cat;

ALTER TABLE contour DROP COLUMN value;

それぞれdumpします。

pgadmin3を使ってテーブルのバックアップをします。

contourテーブルを右クリック-メニューからバックアップをクリック

ダイアログで、ファイル名を指定。contour_dump.sql

PLAINをチェック

PLAINオプションでDB作成をチェックしてOK。

今あるcontourテーブルを削除します。右クリック-削除

suiikiテーブルの場合はファイル名suiiki_dump.sql

データのみをチェック。OK。

出力されたファイルでは、データの挿入先がsuiikiテーブルになっているのでファイルの中でsuiikiとあるところをcontourに置換。

contour_dump.sqlからPostGISに流し入れます。

psql -f contour_dump.sql map

このままでは、gidがプライマリーキーとして設定されていて、suiikiで同じgidを持っているデータが挿入できないので、gidカラムを削除。 キー設定せずにgidというカラムをserial型で作成。データベーステーブルとしてはおかしなものになりますがこのまま使う訳ではないのでいいでしょ う。

(07/09/26 追記 データ構造とGRASSの挙動をあまり理解できていないため、何だか煩雑なことをしていますが、結局は必要なカラムだけを残してマージすればいいみたいです。)

psql -f suiiki_dump.sql map

これで、無事海岸線の入ったcontourが完成。

pgsql2shp map contour -P passwd

カレントディレクトリにcontour.shpが出来ているのでqgisでみてみよう。

• GRASS5で

ベクターとラスターをインポート

加工されたcontour.shpとmosaic.tifが手に入ったので、これからGRASS5での作業に入ります。

GRASS5で適当なlatlonロケーションを作ります。

grass54

ここからはgrassコンソールで作業します。

d.mon x0

v.in.shape input=contour.shp output=contour label=label

ここでは、labelオプションで明示的にラベルをつけるカラムを指定します。

このままではv.digit出来ないのでv.suppurtします(あまり理解できていません)。

v.support map=contour

(07/09/26 追記 ここを省略して、v.build.polylinesが通るならその方がいいと思います。

v.build.polylinesをするとカラム名も新たに作られるみたい。CAT_IDが対象のカラム名。未確認。そうなると、v.surf.rstのzcolumn=ELEVはzcolumn=CAT_IDに変わる。v.infoで確認すること。)

そのままでも、充分きれいにつながった等高線なのですが、欲張ってもっとラインを減らしてみましょう。

g.rename vect=contour,contour_nonpoly

v.build.polylines -c input=contour_nonpoly output=contour

g.remove vect=contour_nonpoly

さて、今度はラスターをインポートしましょう。

r.in.gdal -o input=mosaic.tif output=mosaic

100%

CREATING SUPPORT FILES FOR mosaic.red

SETTING GREY COLOR TABLE FOR mosaic.red (8bit, full range)

100%

CREATING SUPPORT FILES FOR mosaic.green

SETTING GREY COLOR TABLE FOR mosaic.green (8bit, full range)

100%

CREATING SUPPORT FILES FOR mosaic.blue

SETTING GREY COLOR TABLE FOR mosaic.blue (8bit, full range)

あれー、ここでは三つのレイヤーに分解されてしまいました。出力のときにRGBに変換されているんですね。まぁ、ここまでくると裏のデータなんか必要ないので単純に合成してしまいましょう。

r.composite r_map=mosaic.red g_map=mosaic.green b_map=mosaic.blue levels=32 output=mosaic

ここまで来るのに、まる一日はかかってしまいました。数日前におかした間違いを繰り返しています。こうやってメモを残すのは大事ですね。

やっと、v.digitを使う場面がやってきました。

海岸線を1mの標高に指定

v.digitを起動します。

デジタイザーの選択になりますが、接続されていないのでそのままEnter。

編集するファイルの指定は、contourと入力してEnter。

マップの情報は変更する必要がある項目は変更して、＜esc＞＜enter＞。

プロットしていいかと確認が出るのでEnter。

島の海岸線が残っているため広範囲がプロットされます。

目的の地域にズームしますが、

Digitize Edit Label Customize Toolbox Window Help Zoom Quit

とかかれた欄を注目してください。ここがメインのコマンドになります。各コマンドの頭文字を入力すると、そのコマンドのページに移動したり、コマンドが実行されます。

例えば、大文字の「Q」を入力すると、終了していいかと確認が出ます。それに「y」と答えると終了処置がされます。最後に圧縮して保存するかと確認が出ます。それに答えると、v.digitが終了されます。これで、いつでも好きなときに終了できますね。

ここでは、ズームしたいので「Z」を入力します。

Buttons:

Left: Zoom menu

Middle: Pan

Right: Quit

という表示が出ます。ここではマウスのクリックによって動作が決まります。クリックはモニター画面上で行います。

ここで、右クリックするとズームコマンドが終了します。中ボタンをクリックするとクリックされた点を中心に移動します。

左クリックをしてみます。

Buttons:

Left: 1. corner

Middle: Unzoom

Right: Main zoom menu

右クリックで、先ほどのメニュー階層に戻ります。中ボタンでは表示が縮小されます。

左クリックします。

Buttons:

Left: 1. corner (reset)

Middle: 2. corner

Right: Main zoom menu

右クリックで、元の階層に戻ります。左クリックでは、一つ目のコーナーを選び直せます。中ボタンをクリックすると、左クリックで選択した点と、中ボタンで選択した点の範囲が拡大して表示されます。

この、ズームや移動が自由にできると入力作業がはかどります。

さて、拡大してくると表示がよく見えるようになります。

ここで、青く表現されているラインは情報の入っていないものです。

これが、情報を入力されると赤に変化します。

まず海岸線を1mと指定してみましょう。

よく知っている場所に移動して、海岸線を見つけます。

ズームのメニューから抜けます。右上に「Main Menu」と出ている初期画面になっていると思います。ここでコマンド行にある、「Label」を選択します。「L」と入力するんですね。

右上に「Label Menu」と書かれたメニューに変化します。ここでラインに情報を入力する作業をしていきます。

ここでコンソールの真ん中には「Label Menu」で行える作業が列記され、そのコマンドのキーが示されています。

「l - Label Lines」を選んでみましょう。ラインにラベルをつけるコマンドです。

「l」キーを押します。

「Do you wish to enter line labels?[y]」

と聞かれます。ラベルを入力するためにそのコマンドを指定したわけですから、デフォルトの選択どおり「y」です。そのままEnter。

「Enter Category Number (0 to quit):[0]」

入力したい情報を入れます。ここで0を入力するとそのまま終了します。ここでは「1」と入力しましょう。ここで、「0」に指定できないため、便宜上海岸線を1mに指定するわけです。

「Enter a label (＜cr＞ to quit):」

ここではそのラインにつける名前などを入れるんでしょうか。必要ないので＜cr＞(キャリッジリターン=Enter)。

Choose line:

Buttons:

Left: Choose line

Middle: Accept chosen line

Right: Abort/Quit

ここでもマウスを使います。右クリックをすると先ほど指定した情報の入力作業が終了します。海岸線にあたるラインに触れるようにマウスカーソルを移動し左クリックします。ラインが黄色に変化し、正しいラインを選択できたか確認できます。

コンソールの方をみてみましょう。

Line#: 70 Category: (unlabeled)

Line is Not labeled

70番のラインで、これまでラベルされていないと言っています。

ここで、中ボタンをクリックします。赤く変化しました。これで情報が入力されています。

同じラインを左クリックしてみましょう。

Line#: 70 Category: 1

Line is Category 1

70番のラインに1が入力されています。

ここで右クリックをすると、新しい情報の入力に移ります。ここで10を指定し、10mの等高線を指定していくことでも目的は達成できます。

それでは、v.digitの便利な機能を眠らせたままになってしまいますので、別の方法を試します。

ラスターを参照しながら、等高線に標高値を指定

一旦「Label Lines」を抜けましょう。「Enter Category Number (0 to quit):[0]」をそのままEnterします。「Label Menu」に戻ります。

ここで、コンソール右側の中段に「i - Contour interval: < 5>」等高線の間隔が指定されています。「i」キーでこれを変更する事が出来ます。

「i」と入力して質問に「10」と答えます。「i - Contour interval: < 10>」に変化しました。これは、複数の等高線を選択したときにその等高線に指定される数字の間隔は10である、ということです。25000分の1地形図の場合は10m間隔ですので10に変更しました。これは入力作業をする度に必要になりますので、忘れないようにします。

さて、ここからが本番です。「c - Label Contours」複数の等高線にラベルをつけるコマンドです。まず入力してみましょう。「c」

ここでは情報を入力せずに先にラインを選択します。

10mとわかっているラインを選択します。左クリック-中クリック

コンソールに「Enter elevation for this line:」と出ます。「10」と入力します。このラインに10が指定されました。

近くにある30mのラインを選択します。左クリック-中クリック

黄色のラインに青のラインがはさまれた形になります。

「Enter elevation for this line:」に「30」と入力してみましょう。

三本のラインすべてが赤く変化しました。三本をそれぞれ左クリックして情報を見てみましょう。意図したとおりになっていますね。

30mのラインから3本離れた60mのラインを見つけます。そのラインを指定しましょう。

左クリック-中クリック

「60」と入力します。次の指定で先ほど赤く変化した10mのラインを指定してみましょう。情報を入力する前に指定が すんでしまいました。この様にすでに情報を持っているラインは、その情報を他のラインの指定のときに役立てる事が出来ます。赤いラインの間に青いラインが 存在するときは、すでに情報を持っている赤いラインを両側で指定するだけで、キーボードをさわる必要もありません。

ラインの指定をするときに、その経路に同じラインが複数回交差していたり、線をしっかり横切るように指定しないとエラーが出てうまくいかないですが、コンソールに出るエラーを元に対処しましょう。

さて、海岸線から近い等高線は数えることで、その標高がわかりますが、内陸にいくと難しくなってきます。

先ほどインポートしたラスター画像を背景に置くことで、参照しながらの指定作業が出来ます。

「c - Label Contours」の作業を中断しましょう。モニター上で右クリック

「Label Menu」に戻りました。コンソール下部のメニュー項目のうち、「Customize」を選択します。「C」キーを入力します。

「Customize Menu」に入ります。

「B - Select A Backdrop CELL Map None」とあります。背景のラスター画像が「無」になっていますので、指定してあげましょう。「B」を入力。ラスター画像の選択をします。先ほどmosaicがインポートされているのでそれを入力します。

「Do you want to automatically redraw backdrop on re-window? (y/n) [y]」移動してもラスターについてきてほしいので、「y」。

等高線がうまく地図に乗っていますか。先ほど指定した海岸線、等高線は地図上の情報と違っていないですか。実際の作業では、背景がカラーだとライ ンが見にくいので、グレースケールの画像の方がいいでしょう。うまい具合にmosaic.blueなどがありますのでそれを使った方がいいですね。

これで、うまく背景に地形図を出せたので「Customize Menu」から抜けましょう。「q - Return from whence we came」で元いた場所に帰れます。

すでに指定した等高線が間違った標高に指定されているとわかったときは、そのラインのラベルを削除します。

「L - Un-Label Lines」を使います。「L」を入力します。ラインを指定し、確定します。ラインが青く変化しました。

これらの作業を繰り返すことで、すべての等高線に情報を入力します。

ここまでで、疲れたので一旦休みます。

残った作業は、後日詳細を書くかもしれません。

赤い範囲は分析に使えます。定期的にGRASS6にエクスポートした上で、インターポーレートしてGISのデータとして活用してあげましょう。

shape形式でエクスポート

r.out.shape type=line map=contour cats=none prefix=contour.shp

POSTGISに流し込むためにSQLを作成

shp2pgsql -s 4612 -W CP932 -d contour.shp contour>contour.sql

POSTGISに読み込み

psql -f contour.sql map

pgadmin3の中で0とラベルされているラインと、ラベリングされていないラインを削除

DELETE FROM contour WHERE elev=0;

DELETE FROM contour WHERE elev=NULL;

elev=1 の行を0に変更。保存。

(07/09/26 追記

v.surf.rstでは0のデータは無視されるかもしれません。海岸線がおかしいようならelev=-0.001などにしておいたほうがいいかもしれません。)

POSTGISからshapeに書き出し。

pgsql2shp map contour -P passwd

GRASS6にて、インポート。

v.in.ogr -o dsn=contour.shp output=contour

(07/10/20 追記 このままだと、latlonのロケーションにしか読み込めないが、gdalwarpを使ってUTM環境にインポートした方がいいかも。)

g.region 必要な解像度を持ったmap あるいは location/mapset/WIND を直接編集

d.zoom 必要な範囲

v.surf.rst

z値を持つカラムを入力するが、インポートの際に大文字に変換されているので注意。v.info -c contourで確認する。

こんな感じ

v.surf.rst input=contour zcolumn=ELEV elev=elev

mapcalc> demifb0 = if ( elev - 0, elev, 0, 0 )

とすると0m以下の地域が0mで埋められる？

GRASSを用いた地理情報システム入門

07/09/28 追記

石垣島は作成開始から3日で出来上がりました。

この位のクオリティーです。

07/10/17 追記

試行錯誤するには時間のかかってしまうv.surf.rstのテンション、スムーシングのオプションが解説されていました。

Spline Interpolation

FOSS4G06WKSlidar5_helena.pdf

Comments (0)

post 00:10 IST

数値地図 25000(地図画像)をGeoTiffに

2007/09/12

別稿より移動

「数値地図 25000 (地図画像)」をGeoTiffにする

余白の除去

GIMPを使って、地図の周辺余白部分を取り除く。

選択ツールのxをkanri2kw.csvの図郭四隅ﾌｧｲﾙ座標左上Pの数字、yを左上Lの数字、幅を右下P-左上Pの数字、高さを右下L-左上Lの数字とすることで図郭のみが選択される。

選択範囲をコピー、新規画像に張り付け、

図名.tifとして保存、パックビット圧縮。

tfwファイルの作成

yaskeyさんの各種ツールを使わせてもらう。

経度緯度の表示形式変換でkanri2kw.csvの図郭四隅経緯度の左上と右下を変換しておく。

その数値を、Georeference with World file v0.2のそれぞれ、X座標、Y座標に入力。保存しておいたtifファイルの画素数を入力すると計算によりtfwファイルが得られる。図名.tfwとしてtifと同じディレクトリに保存。

grassのr.in.gdalで取り込むことができる。

一図郭につき2ドットほど狂いがでる。台形の画像を長方形にしてしまっているため。

この入力方法をつかっているうちは、「ウォッちず」と解像度にそれほど違いはない。精度は「ウォッちず」の方が上か。「ウォッちず」経由だと図郭の範囲が狭いため同一地域での枚数は多くなる。

オリジナルのtifファイルをそのまま読み込めれば、インデックスカラーで入るので文字の下の等高線情報なども温存できるんでしょうか。

07/09/14 追記

別の方法でうまくいけるようになりました。

create DEM partIII

Comments (0)

post2007,09,09 08:20 IST

compile and install

2007/09/09

GRASSとQGISの最新版を追いかけようとすると、野良ビルドしなくてはいけなくなりますが、aptの恩恵を捨てたくないのでdebuildできるソースを探してきます。

GRASS

cvs -z3 -d:pserver:grass-guest@intevation.de:/home/grass/grassrepository co -P grass6

cd grass6

(07/10/27 追記

cd debian

mv patches/ patches-cvs

svn co svn://svn.debian.org/pkg-grass/packages/grass/branches/6.3/debian .

cd ..

)

debuild -us -uc

cd ..

sudo dpkg -i *deb

QGIS

svn co svn://svn.debian.org/svn/pkg-grass/packages/qgis/trunk

cd trunk

./debian/controlのDependsに含まれるgrass-devをlibgrass-devに変更。(07/10/27 追記 ここは必要なくなったみたい)

debuild -us -uc

cd ..

sudo dpkg -i *deb

unstable環境では、

dpkg-checkbuilddeps: Unmet build dependencies: で要求されたパッケージを入れてあげればビルドできるのかな。こちらでは出来ています。

07/09/22 追記

出来心でqgis 0.9をインストールしてみました。checkinstallを使ってaptitudeでアンインストールできるようにします。

sudo aptitude purge libqgis1 libqgis1-dev qgis qgis-plugin-grass

svn co https://svn.qgis.org/repos/qgis/tags/Preview-0.9.0-1 qgis-0.9.0

cd qgis-0.9.0

mkdir build

cd build

ccmake ..

GDAL includeがずれているので修正。

make

sudo checkinstall

GRASSの中で使ってみたのですが、保存してあるプロジェクトがうまく読めないため削除しました。

また今度やってみよう。

07/10/20 追記

svn co https://svn.qgis.org/repos/qgis/trunk/qgis qgis-0.9.0

cd qgis-0.9.0

./debian/controlに含まれるgrass-devをすべてlibgrass-devに変更。(07/10/27 追記 必要なくなりました。)

debuild -us -uc

dependsには含まれていないですが、python2.5-devが必要な様です。

sudo dpkg -i ../*deb

qgis0.8で保存したプロジェクトはうまく読み込めなかったのですが、新しくレイヤーを読み込んでいくぶんにはうまくいく様です。互換性がないのかな。

これまでのところ、不具合は感じられないのでこのまま使ってみようかと思います。

Comments (0)

post2007,09,05 01:23 IST

create DEM partII

2007/09/05

松田順一郎さんの「LEAVES of GRASS GIS」が404 Not Foundになってしまったので、他人のふんどしでとっていた相撲がとれなくなりました。(07/10/22 追記 今確認したら復活していました。よかった、よかった。)

心を入れ替えて、手順を残しておくことにしました。

記事が完成する時は来るでしょうか。

grass 6.3.0-1

qgis 0.8.1-1

qcad 2.0.5.0-1-2

varicad2007-en 2.01-1

pgadmin3 1.4.3-2.1

野良インストール

GRASS 5.4.1 http://grass.itc.it/grass54/source/grass-5.4.1.tar.gz

varicadはプロプライエタリのソフトですが、試用期間が一ヶ月あったのでその間に使用しました。dwgのファイルをdxfに変換するのに使いました。JGrassのファイル入力メニューにdwgがあるようなのでそれで読み込めるかもしれません。

等高線ベクターの準備

まず、等高線の書き込まれたcadファイルを手に入れます。「数値地図 25000 (地図画像)」などを利用すればジオコーディングされているので手順がだいぶ単純になると思います。

kunigami、oogimi、higashiの三村をまとめてしまいます。が、dxfの中では座標が指定されていないため、xyのロケーションに読み込んだ時点で3つのファイルの位置関係はむちゃくちゃです。

これを座標指定して、隣りあう3つの村にします。

ベクターファイルをラスターに変換します。ベクターのままだとi.pointsが使えないため。

v.to.rast

三つのdxfのうちhigashiでは、そのままのリージョン設定では出来たラスター画像の解像度が低かったため、変換の前にWINDファイルをいじって、解像度を上げました。

これらをlatlonのロケーションにあるベースマップにi.pointsします。

i.group

i.target

i.points

出来たPOINTSファイルを使ってベクターをprojectionします。

v.proj

projectされたベクターをエクスポートします。

v.out.ogr

これを別のロケーションにインポートすればいいのですが、このままでは各ファイルのカラム名がバラバラで簡単にマージできないため、それを統一するために一旦POSTGISに読み込みます。

shp2pgsql

psql

必要なのはz値を入力してあるカラムだけなのでいらないカラムを消去します。

pgadmin3

各テーブルのz値の入ったカラムを同じ名前に変更します。elevなんてどうでしょうか。

各テーブルをダンプします。ダンプされた各sql文のテーブル名を同じ物に変え、droptable文を削除します。

そのsqlを順番に読み込ませれば、POSTGISの一つのテーブルに三村分の等高線情報が入ります。ついでに海岸線のベクターもマージしました。

それをpgsql2shpで出力し、GRASSにインポートすれば広範囲を一度にv.digit出来ます。別々のベクターをv.surf.rstしてラスター画像をマージすると境界や重複部が変な感じになりますがベクターの段階でマージするとマシなのではないでしょうか。

地形図から等高線を抜き出すときはr.patchしてからr.reclass、r.thin、r.to.vectといって、v.out.ogrでgrass5のロケーションに移動すればいいのか。

広い範囲のDEMを作成するとき、皆さんはどうしているんだろう。

GRASS5でv.digit

v.in.ogr

v.support

v.digit

DELETE FROM contour WHERE elev=0

07/09/15 追記

と、ここまで書いてきましたが、数値地図25000(地図画像)を手に入れてみたらそっちの方がいいようなので、その際の手順を別稿で書くことにしました。もうちょっと丁寧にかけるかも。

Comments (0)

post2007,08,20 00:58 IST

orthophoto import and patch

2007/08/20

国土交通省の「オルソ化空中写真ダウンロードシステム」で、ちょっと古い航空写真をオルソ化した画像が公開されています。

「photo import」で扱った画像と同じ物ですが、広い範囲を楽に手に入れられますので使えるようにしてみます。

サイトにアクセスして、必要な地域の画像をダウンロードします。スクリプトなりで自動化できればいいのですが地道に人力でやり遂げました。

google earthに取り込むためのスクリプトが色々と出ているようなので工夫をすればできるんでしょうね。

すべてをunzipして解凍します。svg、jpg、xmlなども含まれていますが、ここでは使いません。

出てきたtifファイルをlatlonのロケーションに読み込みます。

#!/bin/sh

for file in *.tif

do

r.in.gdal -o input=$file output=${file%.tif}

g.region rast=${file%.tif}.red

r.composite red=${file%.tif}.red green=${file%.tif}.green blue=${file%.tif}.blue levels=64 output=${file%.tif}

g.remove rast=${file%.tif}.red,${file%.tif}.green,${file%.tif}.blue

g.remove group=${file%.tif}

done

地域によって違うと思いますが、うちで入手したファイルはすべて「128」から始まっていますので、作られたレイヤーも128からの数字で名前が つけられています。一画像につき一つのレイヤーが作られています。すべてを敷き詰めた絵を見るには手間がかかるので、以下の手順でモザイク画像にします。

MAPS=`g.mlist type=rast sep=, pat="128*"`

g.region rast=$MAPS

r.patch in=$MAPS out=mosaic

07/08/22 追記

すべての画像を一枚にまとめると、表示に時間がかかります。

ハンドリングに支障をきたすためwgsのロケーションに移す際に、四つに分割しました。

r.projで移動する際には、操作時点のregionの範囲に限られるので、regionを設定してから、r.projします。

r.proj input=mosaic location=latlon mapset=hoge output=kunigami1

qgisの中ではpyramidの構築ができるので、作ってしまえば表示は早いです。

と思ったら、保存されてない？

qgisを再度立ち上げると、ピラミッド画像がなくなっていますね。

まぁ、いいか。

オルソ化空中写真ダウンロードシステム

r.patch Manual

Comments (0)

post2006,07,18 15:22 IST

PostGIS

2006/07/18

PostGISも使ってみることにしました。

postgresql-8.1-postgis1 1.1.2-7

postgresql-postgis-common 1.1.2-7

postgresql-postgis-java 1.1.2-7

aptitudeでインストール。

PostgreSQLは既に使っているので、ユーザー登録、パスワードなどの設定はすんでいる。

# mktemplate_gis

postGIS用のテンプレートを作る。

$ createdb.postgis map

SPATIAL_REF_SYSテーブルの入ったmapデータベースができました。

$ shp2pgsql -W CP932 -d kawa.shp river |psql map

で、kawa.shpの内容が入ったriverテーブルができました。

shpはsjis、dbはutf-8ですがちゃんと変換されておりました。

日本語のカラム名はQgisではばけちゃいました。

今持っているベクターデータは全部移した方がいいのかなぁ。

‘07 03 10 追記

マシンが変わったのでPOSTGISの再インストールをしましたが、手順がちょっと変わりました。

$ createdb map

$ /usr/lib/postgresql/8.1/bin/createlang plpgsql map

$ psql -d map -f /usr/share/postgresql-8.1-postgis/lwpostgis.sql

たまに、qgisのオンザフライ投影でうまくいかないことがあるので、shp2pgsqlで-s 4612などのオプションも使います。

$ shp2pgsql -s 4612 -W CP932 -d kawa.shp river> river.sql

$ psql -f river.sql map

できあがり。

Comments (0)

post 15:20 IST

Terrain Analysis

2006/03/14

できあがったbasemapsを使ってどういった分析ができるのか、調べてみます。

勉強に使わせてもらうサイトのリスト。

TERRAIN ANALYSIS AND EROSION MODELING

Terrain modeling and Soil Erosion Simulations for Fort Hood and Fort Polk test areas

Erosion/deposition modeling with USPED using GIS

Landscape soil erosion modeling for spatial conservation planning: GIS-based tutorial

HELENA MITASOVAさんによる例

FREC 682 — Spatial Analysis

デラウェア大学のGISコースシラバス

The GRASS GIS software - GIS Seminar

GRASS-R-Geostatistics

Markus NetelerによるGRASSとRの連携についての解説

“S”に対抗するオープンソース統計解析アプリケーション”R”

GRASS 6 in a nutshell

Slides

Geographic Information Analysis: From GIS to GIA The Science of the System

Kernel Home Range Density

Pachacamac GIS Project: A Practical Application of Geographic Information Systems and Remote Sensing Techniques in Andean Archaeology

The USC GIS Research Laboratory

UC Davis

real distance analysis

RでGIS

環境アセスメントベースマップ整備マニュアル

環境アセスメントベースマップ整備マニュアル

環境省総合環境政策局環境影響評価課

A Spatial Clustering Technique for the Identification of Customizable Ecoregions

GRASS GIS - Simulation Models

The TerraFlow Project

The r.le manual: Quantitative analysis of landscape structures

r.le.setup Manual

r.li is a more flexible and faster replacement of the old r.le.

r.li Manual

FRAGSTATS

FRAGSTATS: spatial pattern analysis program for quantifying landscape structure.McGarigal, K., and B. J. Marks. 1995.

Landscape Analysisという解析

東京大学 高橋さんの「てくてくGIS」内コンテンツ

v.out.kml - A shell script to generate KML files (Google Earth) from GRASS vector layers.

GISによる都市緑地資源の環境価値の評価

JGrass documentation downloads

GRASS GISによる森林ＧＩＳ入門山本一清さんによるチュートリアル データ構造の解析や、Rubyによるデータ変換プログラムの自作など

Power User Workshop: GRASS image processing (with focus on optical and Lidar data) — FOSS4G2006

Spline Interpolation

MultiSpec(c) 概説 日本語訳 マニュアル リモートセンシング解説

Comments (0)

post 15:20 IST

Win and GRASS

2006/03/12

Windows版のGRASSとQgisが出ているようです。

GRASS for Windows

使ってはいないのですが、日本語が通ればQgisをインターフェイスとしたGRASSの利用が広まってくれるかもしれないですね。

bittorentoを利用してダウンロードするように変更されましたね。

(06/06/30 追記 vmware上のWin98で使ってみました。日本語通りますね。素晴らしいですね。Win上でのGRASSの使いかたがわからない。もうちょっと触ってみます。)

(07/08/22 追記 使えるようになっていません。Windowsを立ち上げることがほとんど無いので、触る機会がありません。)

(07/09/27 追記 この通り実行すればいいみたいです。grass Win)

この2ヶ月ほどで、GRASSを使う上でのGUIの改善がすごい勢いで進んでいるように思います。

以下、現状のGRASS cvs とQgis Ver.0.8 のスクリーンショットです。

スタートアップ画面。

新しいロケーションの作成時、EPSGコードで指定して投影法などを選ぶことができるようです。

Qgisの日本語対応状況。

以前に較べると、メニューの日本語が少なくなったようですが、しばらくすると対応されてくると思います。表示など、全然問題なく行えます。

画面は、shpファイル中のsjisの地名を表示させています。

grass-toolsのmapbrowserです。表示の選択や、ちょっとした操作もできます。

grass-toolsのmodulesです。

かなりのmoduleが操作可能になっていてびっくりしました。細かいコマンドラインオプションの指定はできませんが、気軽に使えるので便利です。

僕は、商用のソフトを使ったことがないので、それらと比べてどの程度の差があるのかわかりません。しかし、これらのバージョンを元にした日本語での解説が増えれば、日本での使用者が劇的に増えるんじゃないかと思います。

07/08/22 追記

QuantumGIS User Guide 日本語ユーザーガイド

Comments (0)

post 15:19 IST

DEM Export

2006/03/12

DEMをエクスポートしてみました。

AAIGrid Int16の形式です。

Arcなどで無事に読めるでしょうか。

2007/08/20 追記

Arcを持っていないため検証できていませんが、以下のサイトに答えが書いてありました。

Exporting GRASS Raster Data to ArcGIS CA Soil Resource Lab

2007/09/03 追記

ラスター画像をgrass5に持っていくと色がおかしいので、以下のコマンドでだしました。

r.out.tiff input=map output=map.tif -t compression=deflate

これでDEMとして使えるのだろうか。

人に渡す必要がある時は、上のサイトの情報でやってみよう。

Comments (0)

post 15:18 IST

photo import

2006/02/25

国土画像閲覧システムから、目的の画像をダウンロードし、若干加工してxy座標系に読み込みます。

その後、同地域をカバーするDEMを使ってオルソ幾何補正して、目的のロケーションに読み込みます。

使うコマンド

r.in.gdal

r.composite

g.region

i.group

i.target

i.ortho.photo

i.ortho.photo currently requires the elevation model to be in meters, and the target location to be a standard (e.g. UTM) coordinate system.

地形にあわせて変形された結果、長方形だった画像もこういう感じにエッジまで変形されます。

nvizで、DEMに乗っけるとこんな感じに見えます。

ESRI オルソ幾何補正とは？

オルソフォト

単画像オルソ幾何補正

ERDAS IMAGINE オルソ幾何補正機能の処理時間と品質の検証

単画像オルソ幾何補正 プレゼンテーション

HotTopic3： IMAGINE OrthoBASEによる米軍撮影航空写真のオルソ幾何補正

i.ortho.photo Manual 機械翻訳

Comments (0)

post 15:17 IST

create DEM

2006/01/12

デジタル標高モデルを自力で作ります。

参考にした松田順一郎さんの「LEAVES of GRASS GIS」そのままなので、そちらも参照してください。

grass5を起動。

使うコマンド

r.in.gdal

r.reclass

r.thin

r.line

v.support

v.digit

v.surf.rst

i.group

i.target

i.points

r.proj

i.rectify

目的の地形図を入手しインポートします。

http://watchizu.gsi.go.jp/

4隅の緯度経度情報も後ほど必要になります。

ここでできあがったDEMは、utmのロケーションでr.mapcalcを使って接合します。

現在、この程度できています。地形図4枚分です。

結構大変な作業で、時間がかかります。

Stereo PhotogrammetryでDEMの自動生成ができれば楽なのですが、現時点ではオープンソースでのその用途のソフトはうまく動かせていません。

stereo-grass、e-foto e-fotoに期待をしたいところです。

いろいろ試してみました。

国頭全域の等高線情報の入った、DXFファイルを手に入れました。

かなり大きいのですがcvsのGRASSでは読みこめたので、これを元にDEMを作成する作業をしました。

まず、xyのロケーションに読みこみます。

v.in.dxf

Number of lines : 171961

こんなにラインがあってはたまらないので、v.cleanしてみます。

v.clean

Number of lines : 66140

かなり減らすことができました。

短いラインを消す作業をしたので、崖の表示などがさっぱりしました。赤いのが消されたラインです。

これを、utmのロケーションにいったん読みこみます。

v.transform

v.out.ogr

v.in.ogr

utmで海岸線のベクターデータとマージします。

v.patch

v.out.ogr

grass5でv.digitして上記と同じ作業をします。

等高線が癒着したり、細切れになったりしていないので作業は楽になります。

作業開始から数時間で、赤色の範囲が処理できました。

1/25000地形図からDEMを作る．

(07/09/04 追記 リンク切れしてしまいました。)

ERDAS IMAGINE オルソ幾何補正機能の処理時間と品質の検証

07/08/30 追記

Trace vector contours from a scanned map

07/09/03 追記

v.digitのバックグラウンドをラスター画像にする方法がわかりました。

C(Customize)-B(Select A Backdrop CELL Map)

で、画像選択になります。

座標指定をしてあって、ちゃんと重なる画像があれば下のように効率的に標高指定が可能です。

いやー、楽だ。

Comments (0)

post 15:15 IST

import BASE MAP

2005/12/29

windows版景観閲覧用フリーソフトであるカシミールを使って、広範囲の接合された地形図を作り、xyロケーションに読み込みます。

その後、目的のロケーションに読み込みます。

使うコマンド

r.in.gdal

i.group

i.target

i.rectify

カシミールで地形図を作る

カシミールの起動画面で出てくるダイアログの中の『「ウォッちず」を利用する。』のボタンをクリックすると国土地理院のサーバーから地形図が読み込まれ、自動的に正規化、接合が行われます。

目的の範囲が全て読み込まれたら「編集」-「選択範囲を決める」。地図上をドラッグすると選択されます。

ドラッグ開始地点の緯度経度、ドラッグ終了地点の緯度経度を記録しておきます。

あるいは、「編集」-「選択範囲を緯度経度で指定する」を選んでもいいです。

「ファイル」-「表示画像を保存」-「選択範囲を保存」します。

このとき「表示測地系」はwgs84を選択していたと思います。

この画像は「ウォッちず」 から直接ダウンロードした画像に比べると劣化しているようですので、DEMの作成には適しません。DEMを作る場合は地形図一枚ずつを別に読み込むように しています。あくまで、他の画像を読み込む際のGCP(グランドコントロールポイント)を設置するためのものと考えています。

画像の読み込み

事前にGIMPで、8色インデックスカラーに変換してもいいですし、pngのまま読み込んで3色に分解されたものをr.compositeしてもいいです。

xyのロケーションにr.in.gdalを使って読み込みます。

r.in.gdal input=~/doc/gis/data/mapname.png output=mapname

g.list rast で確認すると mapname の出来たことがわかります。このとき、mapname.blue mapname.green mapname.red が出来ている場合は、3チャンネルに分解されているのでr.compositeする必要があります。

g.region rast=mapneme

d.rast mapneme

で、図幅全体をみることが出来ます。d.zoomで拡大、縮小が可能です。

全景

座標の指定

地形図の場合、きちんと正規化されているので4隅の座標がわかっていればその内側全てが正しい位置関係になっています。写真などの場合は1枚の図の中でも細かい歪みがあるので、地形データを元にして、図を細かく変形しないときちんと地図に乗りません(オルソ化)。

他のロケーションにデータを移動するときには、グループ単位で行います。1枚の画像を移動するときでも、グループを作りそこに移動する画像をいれます。

i.group group=groupname input=mapname

そのグループを移動する目標のロケーションを指定します。

i.target group=groupname location=target mapset=username

特に指定しなければusername=マップセット名になっていると思います。

まず、latlonに読み込ませます。

カシミールで切り出したときの緯度経度情報のフォーマットが、DDD度MM分SS.S秒で表現されているのでDDD.DDDDDDDDDの形式に直します。

緯度経度情報をMANDARA用に変換する方法を参考に4隅の数字を記録します。

yaskeyさんのwebスクリプトが便利ですね。

経度緯度の表示形式変換

準備した画像では

左下 北緯26度42分12秒, 東経128度08分00秒

26.703333 128.133333

左上 北緯26度52分50秒，東経128度08分00秒

26.880556 128.133333

右下 北緯26度42分12秒，東経128度20分30秒

26.703333 128.341667

右上 北緯26度52分50秒, 東経128度20分30秒

26.880556 128.341667

という結果になりました。

これを元にPOINTSファイルを作ります。

$GISBASE(例えば/home/user/gisdata/grass5)/mapset名/group/groupname/POINTS

が、このグループの移動先での座標を指定するPOINTSファイルです。

内容は以下の様になります。

# image target status

# east north east north (1=ok)

#

-5333 -5443 128.133333 26.703333 1

-5333 0 128.133333 26.880556 1

0 -5443 128.341667 26.703333 1

0 0 128.341667 26.880556 1

表のレイアウトの関係で空白が全角文字になっていますが、半角かタブにしてください。

各行の1番目の数字、2番目の数字は元になる図のx座標、y座標です。(web上の情報と違って右上が原点になるようです。Ver.が上がって変わったのかなぁ。)例に出している図の大きさが、5333x5443の画素数がある絵なのでこんな感じになります。

行の最後にある「1」はこの行をデータとして使用する、という意味のようです。

準備ができたら、i.rectifyで目的のロケーションに送り込みます。

i.rectify group=this_group extension=map order=1

完了したらメイルが届くという表示がされますが、届いたり、届かなかったりします。

目的のロケーションで、g.list rastしてグループ内のファイル名+map(extensionで指定した文字列)が存在していれば完了です。

latlonにちゃんと入ったら、今度はutmに読み込みます。

r.proj input=mapname+map location=latlon mapset=username output=mapname

何かshapeファイルで提供されているベクターデーターと重ねてみて、ズレが無いかどうか確認してみましょう。

d.mon x0

g.region rast=mapname

d.rast mapname

d.vect road color=red

ばっちりでした。

LEAVES of GRASS GIS

カシミール３Ｄ

緯度経度情報の変換

Comments (0)

post 15:14 IST

location

2005/12/28

GRASSは、一つのロケーションでは事前に定めた一つの測地座標系しか扱えません。

なので、種々の測地座標系を設定したロケーションを複数作って、使い分けます。

僕は、

xy

grass5

latlon

utm

とそれぞれ名前を付けたロケーションを使っています。

xy

座標の指定されていない画像データをGRASSに取り入れるために使います。ここでlatlonの座標を指定して、i.rectifyで貼り込みます。

grass5

Ver.6のv.digitはうまく使いきれなかったので、Ver.5を使います。その際、サポートしていないprojなどぐちゃぐちゃになるか と思い、特別にロケーションを用意しました。結局、i.rectifyするときなどVer.6でも触っていますが大丈夫なので、xyでやってしまってもい いかもしれません。

latlon

環境省や国土交通省、国土地理院などから発表されている各種データが、平面直角座標系のshapeファイルなので、それらを取り込むときに使います。

最初はellipsoidをBesselと指定していましたが、JGD2000のデータを読み込むと数メートルのズレがでるため、ロケーションを作りなおしました。

datumの指定は無し、ellipsoidはgrs80だと、きれいにおさまります。

こっちをメインに使おうと思っていたのですが、nvizでの表示でバグがでるのでデータのinportで中継に使うだけになってしまいました。

NVIZ -TODO and BUGS

“nviz does not work well with latitude-longitude projections because it

uses decimal degrees instead of metres for the resolution. Thus is why

the vertical exaggeration is such a small number.”

cat latlon/PERMANENT/DEFAULT_WIND

proj: 3

zone: 0

north: 27:53:09.472N

south: 24:02:44.217N

east: 131:19:55.608E

west: 122:56:01.151E

cols: 2141687

rows: 1182010

e-w resol: 0:00:00.014117

n-s resol: 0:00:00.011696

top: 1

bottom: 0

cols3: 2141706

rows3: 1182050

depths: 1

e-w resol3: 0:00:00.014117

n-s resol3: 0:00:00.011696

t-b resol: 1

g.proj -p

-PROJ_INFO————————————————-

name : Latitude-Longitude

proj : ll

ellps : grs80

-PROJ_UNITS————————————————

unit : degree

units : degrees

meters : 1.0

utm

メインの作業スペースになります。germinのGPSなんかだと、データはここに直接入るのかな。

一旦latlonに取り入れた画像はr.projの自動変換機能を使って取り込みます。

ここだとnvizもきれいに表示できます。

cat utm/PERMANENT/DEFAULT_WIND

proj: 1

zone: 52

north: 1

south: 0

east: 1

west: 0

cols: 1

rows: 1

e-w resol: 1

n-s resol: 1

top: 1

bottom: 0

cols3: 1

rows3: 1

depths: 1

e-w resol3: 1

n-s resol3: 1

t-b resol: 1

g.proj -p

-PROJ_INFO————————————————-

name : Universe Transverse Mercator

datum : wgs84

towgs84 : 0.000,0.000,0.000

proj : utm

ellps : wgs84

a : 6378137.0000000000

es : 0.0066943800

f : 298.2572235630

zone : 52

-PROJ_UNITS————————————————

unit : meter

units : meters

meters : 1.0

Grassのデータベース座標系の設定

Bubble://別館/ - B-Wiki - SRIDList

Comments (0)