[OSGeoJapan-discuss] QGISのpythonコンソールで面積を計算してその結果を列に加える方法

2015年 5月 20日 (水) 20:06:39 PDT

飯島さま、

蛇足ながら…もっともローテクな方法を。

ポリゴンの面積計算はベクタ→ジオメトリツール→ジオメトリカラムの出力/追加で
計算できます（ポリゴンの場合、面積と周長のカラムが追加されます）。

プロセッシングツールでバッチ処理をすることも可能です。

Processing→Toolbox→開いたウィンドウ下部で[advanced interface]を選択
→QGIS geoalgrithms→Vector table tools
→Export/Add geometry columnsを右クリックし[バッチプロセス…]を選択
→Input layer, output layerに入力
（Calculate usingは面積計算に使われる測地系です。レイヤの測地系が
緯度経度の時はプロジェクトの測地系をUTMや直角座標系に設定し、
[ProjectCRS]を選択したりしてます）

プログラムが苦手な方に処理をお願いする際など便利かもしれません。。

久保まゆみ


2015年5月21日 11:47 Hayato IIJIMA <hayato.iijima ＠ gmail.com>:

> 垂水様、皆様
>
> 山梨県の飯島です。次々と方法をご提示いただき、ありがとうございます！さすがOSGeoメーリングリスト......
>
> 取り急ぎ、今木様の方法でやってみたらできましたので、ご報告させていただきます。
>
>
> 先ほどのsample.shpがある環境でRを起動して、以下のように実行しました。あといまさらですが、最初に添付したsample.shpのEPSGコードは2450です（すいません、大事な情報を忘れていました）。
>
> setwd("C:...") #SHPがある場所へのパス
> library(maptools)
> crs <- CRS("+proj=tmerc +lat_0=36 +lon_0=138.5 +k=0.9999 +x_0=0 +y_0=0
> +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
> hrr.shp <- readShapePoly("sample", verbose=TRUE, proj4string=crs)
> library(rgdal)
> hrr.shp.2 <- spTransform(hrr.shp, CRS("+init=epsg:2450"))
> library(rgeos)
> AREA <- gArea(hrr.shp.2, byid=TRUE)
> library(foreign)
> d <- read.dbf("sample.dbf")
> d$AREA <- AREA
> write.dbf(d, file="sample.dbf")
>
> これを後は複数ファイルで処理できるようにforでループ化すればいいだけです。
>
> 取り急ぎ、お礼まで。
>
>
>
> 2015年5月21日 11:35 株式会社カーネル 垂水 <tarumi ＠ kk-kernel.co.jp>:
> > 垂水です
> >
> > Shapeファイルそのままではないですが、
> > spatialiteにインポートしてSQLで計算もありではないでしょうか？
> >
> >
> > On 2015/05/21 11:31, Taro Matsuzawa wrote:
> >> 松澤です。
> >>
> >> この手のものはもしPythonを使うならQGISよりも
> >> ogrで直接処理したほうがよいのかなと思います。
> >> ipythonとか使えばデバッグも楽ですし、
> >> 何よりQGISのビルドがコケて環境破壊した時でも使えます。
> >> # MacPortsあるあるなのでWindowsならそんな事起きないと思うけど
> >>
> >> 1. 面積計算
> >>
> https://pcjericks.github.io/py-gdalogr-cookbook/geometry.html#calculate-the-area-of-a-geometry
> >> 2. Shapefileに属性テーブルを追加してデータを投入する例
> >>
> https://pcjericks.github.io/py-gdalogr-cookbook/vector_layers.html#create-a-new-shapefile-and-add-data
> >> (新規でファイル作ってるけど、やってることは同じだと思う)
> >>
> >> これらの組み合わせでいける気がします。
> >> あと、Pythonでは他にもいろいろライブラリがあるみたいなので、
> >> 探してみるとよいかもしれません。
> >>
> >> On 2015/05/21 11:12, Hayato IIJIMA wrote:
> >>> 今木様、皆様
> >>>
> >>> 山梨県の飯島です。早速のお返事ありがとうございました。
> >>>
> >>> こんなパッケージがあるとは知りませんでした。Rでできるならその方が得意なので、ありがたいです。早速試してみます。
> >>>
> >>> 他の皆様も別な技をご存知でしたら、よろしくお願いします。
> >>>
> >>>
> >>> 2015年5月21日 10:45 Hiroo Imaki <hiroo ＠ angeli.org>:
> >>>> 飯島様、
> >>>> お久しぶりです。Rでやっても良いのではないでしょうか？
> http://www.inside-r.org/packages/cran/rgeos/docs/gArea
> が参考になると思います。spパッケージを使ってオブジェクトを作る必要がありますが。
> >>>> 参考までに。
> >>>> いまき
> >>>>
> >>>> 2015-05-20 19:48 GMT-04:00 Hayato IIJIMA <hayato.iijima ＠ gmail.com>:
> >>>>> OSGeoメーリングリストの皆様
> >>>>>
> >>>>> 山梨県の飯島と申します。このたびはお世話になります。
> >>>>>
> >>>>>
> >>>>>
> QGISのpythonコンソールで、SHPファイル（添付したようなファイル）のポリゴンごとの面積を計算する方法をご存知でしたら、ご教示いただければと思います。
> >>>>>
> >>>>>
> >>>>>
> ファイルが1つであればSHPファイルを開き、「属性テーブルを開く」→「フィールド計算機」で新しい列を作りその値として計算した面積を入れればいい話です。しかし、このようなファイルが多数ある場合、GUIによる作業は時間がかかるため、pythonで一気に処理できないかと思っています。
> >>>>>
> >>>>>
> >>>>>
> ネット上でQGISのpythonでの処理方法をいろいろ調べてみて、面積を計算して単に表示するとか、新しい列を作るといったことはできるのですが、面積を計算してその結果を新しい列に加えるということがどうしてもできないのです。
> >>>>>
> >>>>> from PyQt4.QtCore import *
> >>>>> from qgis.analysis import *
> >>>>> input_dir="C:/..." #作業するSHPファイルまでのパス
> >>>>> import os
> >>>>> import glob
> >>>>> import processing
> >>>>> files = glob.glob(os.path.join(input_dir, r'*.shp'))
> >>>>>
> >>>>> #面積を「表示」
> >>>>> for i in range(0, len(files)):
> >>>>> vlayer = QgsVectorLayer(files[i], "temp", "ogr")
> >>>>> #Gmail上でタブによる字下げがうまくいきませんが、実際にはタブで行頭を字下げしています（以下同様）。
> >>>>> for f in vlayer.getFeatures():
> >>>>> print f.geometry().area()
> >>>>>
> >>>>> #新しい列を加える
> >>>>> #カラムを加える（試行錯誤）
> >>>>> for i in range(0, len(files)):
> >>>>> vlayer = QgsVectorLayer(files[i], "temp", "ogr")
> >>>>> #Gmail上でタブによる字下げがうまくいきませんが、実際にはタブで行頭を字下げしています（以下同様）。
> >>>>> caps = vlayer.dataProvider().capabilities()
> >>>>> if caps & QgsVectorDataProvider.AddFeatures:
> >>>>> res = vlayer.dataProvider().addAttributes([QgsField("AREA",
> >>>>> QVariant.Double)])
> >>>>>
> >>>>> これらを合わせるとうまくいきそうな気もしますが、今のところできていません。
> >>>>> ちなみに、作業環境は以下のようです。
> >>>>>
> >>>>> OS: Win7 64 bit
> >>>>> QGIS: 2.8
> >>>>>
> >>>>> 調査に出ているためご教示いただいた後の返事が遅くなるかもしれませんが、よろしくお願い致します。
> >>>>>
> >>>>>
> >>>>> --
> >>>>> **********************************************************
> >>>>> 山梨県森林総合研究所 森林研究部 森林保護科
> >>>>>
> >>>>> 研究員 飯島 勇人（いいじま はやと）
> >>>>>
> >>>>> 〒400-0502
> >>>>> 山梨県南巨摩郡富士川町最勝寺2290-1
> >>>>> TEL:0556-22-8005（直）8001（代）
> >>>>> HP: http://www.pref.yamanashi.jp/shinsouken/
> >>>>> E-mail: hayato.iijima ＠ gmail.com
> >>>>> **********************************************************
> >>>>>
>
-------------- next part --------------
HTML$B$NE:IU%U%!%$%k$rJ]4I$7$^$7$?(B...
URL: <http://lists.osgeo.org/pipermail/osgeojapan-discuss/attachments/20150521/98718a7a/attachment-0001.html>