前回の続き
http://qiita.com/motokazu/items/aaa5970ab34044237cc7
前回は特定の地点からの半径で、円を描いてその中に点を生成したが、ある地域だけに生成したいなどなった場合には、円ではだめ。
ということで、国土数値情報のShapeファイルを使って東京内に点を生成します。
ダウンロードページを開いたら、JPGISを選択して、「行政区域」->「東京」を選択、現状で最新の平成27年を選択。アンケートに記入してダウンロードする。
shapeファイルをQGISで開いてみる
東京が表示された！
このshapeファイルを使って、東京に点を打つのである。
処理イメージとしては、
1. 前回のツールで点を打つ
2. shapeファイルを読み込んでおいて shape ファイル内にある点かをチェック
3. shapeファイル内にある点ならOKとする
shapeファイルを読み込んで、使えるようにする
こちらを参考に。
http://stackoverflow.com/questions/7861196/check-if-a-geopoint-with-latitude-and-longitude-is-within-a-shapefile
まずはインストールして使えるように
Shapeファイルの読み込みには fionaを使う
$ pip install shapely
$ pip install fiona
ライブラリを読み込み
# shape
import fiona
from shapely.geometry import Point,asShape
shapeファイルを取り出して、Latitude, Longitudeが示す点が shape ファイルにあるかを確認
with fiona.open(shapefile) as fiona_collection:
shapefile_record = fiona_collection.next()
shape = asShape(shapefile_record['geometry'])
geolat = fake.geo_coordinate(center=centerlat , radius=latradius)
geolong = fake.geo_coordinate(center=centerlong, radius=longradius)
point = Point(geolong, geolat)
if shape.contains(point):
全くcontaints がTrueにならない... と思ったら、shapeファイルが複数のpolygonで構成されている上に、fiona_collection.next() で取ってきているのが1つの区画だけだからだ...きっと
こちらを参考にして、複数のpolygonを融合。助かります！
https://sites.google.com/site/qgisnoiriguchi/vector01/09
QGISを使って結合します。
メニュー「ベクタ」->「空間演算ツール」 -> 「融合」を選択
生成するshapeを選択して、、OK
一つなった！
このshapeファイルを使って点を生成してみる
明らかにcontaintsがTrueになる数が変化した。
最終的なコードはこちら
https://gist.github.com/motokazu/a1cb634aa0d6726039bd
使い方例
$ python genpointsbygeo.my.py --shapefile Tokyo.shp --samples 100 --radius 50000