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

2015年 5月 20日 (水) 19:47:17 PDT

垂水様、皆様

山梨県の飯島です。次々と方法をご提示いただき、ありがとうございます！さすが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
>>>>> **********************************************************
>>>>>
>>>>> **********************************************************
>>>>> Hayato Iijima (Ph.D.)
>>>>>
>>>>> Yamanashi Forest Research Institute
>>>>>
>>>>> 2290-1, Saishoji, Fujikawa, Yamanashi 400-0502, Japan
>>>>> TEL:0556-22-8001
>>>>> HP: http://www.pref.yamanashi.jp/shinsouken/
>>>>> E-mail: hayato.iijima ＠ gmail.com
>>>>> **********************************************************
>>>>>
>>>>> _______________________________________________
>>>>> OSGeoJapan-discuss mailing list
>>>>> OSGeoJapan-discuss ＠ lists.osgeo.org
>>>>> http://lists.osgeo.org/mailman/listinfo/osgeojapan-discuss
>>>>
>>>>
>>>>
>>>> --
>>>> Hiroo Imaki
>>>> Pacific Spatial Solutions, LLC
>>>> 1523 Chatham Colony Ct.
>>>> Reston, VA  20190
>>>> hiroo.imaki ＠ pacificspatial.com
>>>> http://www.pacificspatial.com
>>>> http://www.geopacific.org (GIS info site)
>>>> hiroo ＠ angeli.org (private)
>>>>
>>>
>>>
>>
>
> _______________________________________________
> OSGeoJapan-discuss mailing list
> OSGeoJapan-discuss ＠ lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/osgeojapan-discuss

-- 
**********************************************************
山梨県森林総合研究所　森林研究部　森林保護科

研究員　　飯島　勇人（いいじま　はやと）

〒400-0502
山梨県南巨摩郡富士川町最勝寺2290-1
TEL:0556-22-8005（直）8001（代）
HP: http://www.pref.yamanashi.jp/shinsouken/
E-mail: hayato.iijima ＠ gmail.com
**********************************************************

**********************************************************
Hayato Iijima (Ph.D.)

Yamanashi Forest Research Institute

2290-1, Saishoji, Fujikawa, Yamanashi 400-0502, Japan
TEL:0556-22-8001
HP: http://www.pref.yamanashi.jp/shinsouken/
E-mail: hayato.iijima ＠ gmail.com
**********************************************************