このたびArcGISコミュニティを発足することになりました。
ArcGISは便利なエンジンですが、世に出回る情報量の少なさが玉に瑕ですよね（あったとしても英語ばかりといったことも多々）。ですので、困った時に質問や情報収集が出来るコミュニティの存在が必要だと感じていました。このFacebookグループが日本におけるArcGISコミュニティの先駆けになれればと思っています。
■対象者
ArcGISに興味のある方全員
■対象技術
・ArcMap
・ArcObjects
・ArcPy
・その他
ArcGISに関わるすべての技術を対象とします。またトピック的なものも含めて情報交換できればと思っています。
■今後の展望
・定期交流会の開催（業種を問わず情報交換できればと思っています。学生さんも大歓迎です）
・ArcGISハッカソンの開催
ブログ主はまだまだArcGISを勉強中の身ですので、もっとArcGISを覚えたいという方と一緒に技術を高められれば幸いです。お気軽にご参加ください。
さて、本日はpandasでHTMLの表をスクレイピングしてみようと思います。スクレイピングだったらBeautifulSoupやlxmlなどが有名ですが、HTMLの表をスクレイピングするのは割と面倒くさい作業です。「気象庁の震度データベース」というちょうどいいサンプルがあったので、今回はpandasを使ってみようと思います。←これがとっても便利なんです
ここで過去10年間で起きた震度5強以上の地震を検索してみたいと思います。
52件もあることがわかりました。やはり日本は地震大国ですね。そして、地震の検索結果が表になっていることがわかります。ここから緯度経度を取得して地震のポイントをマップにプロットしようと思います。
このページでもポイントをプロットしています。地震の深さでポイントの大きさを変えていますね。
このページをいったんHTMLとしてローカルに保存します（動的に作られているページなので）。そのHTMLの表をスクレイピングしてみます。
pandasだけでなくlxml、html5lib、BeautifulSoup4もインストールする必要があります。
HTMLの表をスクレイピングして結果をCSVに出力するサンプルです。
# -*- coding: utf-8 -*- import pandas as pd #保存したHTML html = 'D:\python\earthquake\earthquake.html' #HTMLを読込 dataframes = pd.io.html.read_html(html) #表の部分を取得 table = pd.DataFrame(dataframes[5]) table.columns = ['No','Time','Name','Latitude','Longitude','Depth','Magnitude','SeismicIntensity'] #CSV出力 table.to_csv('D:\python\earthquake\earthquake.csv' ,encoding = 'Shift-JIS')
まず、以前と同じようにESRIから全国地図をダウンロードします。
全国市区町村界データ | 製品 | ESRIジャパン
ArcMapで読み込みます。
ポイントのフィーチャークラスを作成します。
フィーチャークラスの属性を以下の地理座標系にします。
「サンプルコード①」で出力したCSVを元にマップ上にポイントをプロットするサンプルです。色々やっているので、要点だけ書いておきます。
1.CSVから読み込んだ緯度経度を60進数から10進数に変換
2.「1」の緯度経度からポイントをプロット
3.「earthquake」フィーチャクラスにカラム追加(Arcpy)
4.「3」のカラムにCSVから読み込んだ情報を付与(Arcpy)
# -*- coding: utf-8 -*- import arcpy import pandas as pd pointGeometryList = [] def main(): #CSV読込 table = pd.read_csv(r"D:\python\earthquake\earthquake.csv",encoding="SHIFT-JIS") data_list = [] for key,rowD in table.iterrows(): if key > 0 : #緯度経度を度分秒に分解 latDo = rowD.Latitude[0:rowD.Latitude.find(u"°")] lonDo = rowD.Longitude[0:rowD.Longitude.find(u"°")] latFun = rowD.Latitude[rowD.Latitude.find(u"°")+1:rowD.Latitude.find(u".")] lonFun = rowD.Longitude[rowD.Longitude.find(u"°")+1:rowD.Longitude.find(u".")] latByou = rowD.Latitude[rowD.Latitude.find(u".")+1:rowD.Latitude.find(u"′")] lonByou = rowD.Longitude[rowD.Longitude.find(u".")+1:rowD.Longitude.find(u"′")] #緯度経度を60進数から10進数に変換 lat,lon = latLong(latDo,lonDo,latFun,lonFun,latByou,lonByou) data_list.append([rowD.No,rowD.Time,rowD.Name,lat,lon,rowD.Depth,rowD.Magnitude,rowD.SeismicIntensity]) CreatePont(lat,lon) table2 = pd.DataFrame(data_list) table2.columns = ["No","Time","Name","Latitude","Longitude","Depth","Magnitude","SeismicIntensity"] arcpy.CopyFeatures_management(pointGeometryList, r"C:\ArcPySample\test.gdb\earthquake") arcpy.env.workspace = ur"C:\ArcPySample\test.gdb" #earthquakeレイヤにカラムを追加 arcpy.AddField_management("earthquake", "Time", "Text") arcpy.AddField_management("earthquake", "Name", "Text") arcpy.AddField_management("earthquake", "lat", "Double") arcpy.AddField_management("earthquake", "lon", "Double") arcpy.AddField_management("earthquake", "Depth", "Text") arcpy.AddField_management("earthquake", "Magnitude", "Text") arcpy.AddField_management("earthquake", "SIntensity", "Text") for key,rowD in table2.iterrows(): cursorJ = arcpy.UpdateCursor("earthquake") for rowJ in cursorJ: if rowD.No == rowJ.OBJECTID: rowJ.setValue("Time", rowD.Time) rowJ.setValue("Name", rowD.Name) rowJ.setValue("lat", rowD.Latitude) rowJ.setValue("lon", rowD.Longitude) rowJ.setValue("Depth", rowD.Depth) rowJ.setValue("Magnitude", rowD.Magnitude) rowJ.setValue("SIntensity", rowD.SeismicIntensity) cursorJ.updateRow(rowJ) #ポイントを配置 def CreatePont(lat,lon): global pointGeometryList point = arcpy.Point() point.X = lon point.Y = lat #日本地図の空間参照を利用 dataset = r"D:\python\earthquake\japan_ver80.shp" spatial_reference=arcpy.Describe(dataset).spatialReference pointGeometry = arcpy.PointGeometry(point,spatial_reference) pointGeometryList.append(pointGeometry) #緯度経度を60進数から10進数に変換 def latLong(latDo,lonDo,latFun,lonFun,latByou,lonByou): lat1 = float(latDo) lat2 = float(latFun) / 60.0 if type(latByou) == float: lat3 = float(latByou) / 60.0 / 60.0 else: lat3 = float(0) lat1 = lat1 + lat2 + lat3 long1 = float(lonDo) long2 = float(lonFun) / 60.0 if type(lonByou) == float: long3 = float(lonByou) / 60.0 / 60.0 else: long3 = float(0) long1 = long1 + long2 + long3 return lat1,long1 if __name__ == '__main__': main()
ポイントがプロットされました。位置もあっているっぽいですね。
ポイントに属性も入っていますね。
次にシンボルの設定をします。個別値として設定します。←ここもArcpyでできるようにESRIさんお願いします。
震度でポイントの色分けを行うサンプルです。
# -*- coding: utf-8 -*- import arcpy mxd = arcpy.mapping.MapDocument("current") lyr = arcpy.mapping.ListLayers(mxd, "earthquake")[0] if lyr.symbologyType == "UNIQUE_VALUES": lyr.symbology.valueField = "SIntensity" lyr.symbology.addAllValues() arcpy.RefreshActiveView() arcpy.RefreshTOC() del mxd
ポイントの色分けができました。次にシンボルの大きさを分けようと思ったのですが、どうやらこれもArcpyではできない模様。うー貧弱すぎる・・・
とりあえずArcMapで変更してみました。いい感じになりました。このようにデータを引っこ抜けば色々な角度から分析できるので便利ですね！
意外なことに都心部では過去10年間に震度5強以上の地震は起きてないようですね。また、東日本大震災の震源はけっこう陸から離れているのにあの津波が起きたわけですから、規格外の大きさの地震ということがわかります。
結果を見ると全国的に大きな地震は起きているのですが、西日本は少ないことがわかりました。沖縄に関しては0です。この結果にプレートとの関係性を絡めたらもっと面白い結果が得られそうです。
とりあえず、今回はここまでとして、次回以降はさらに深く掘り下げた分析をしてみたいと思います。
以上、本日は「 【GISとPython】pandasで気象庁の震度データベースをスクレイピングして震源をマップにプロットする」でした。
さて、本日は神戸市のオープンデータを利用してみようと思います。阪神・淡路大震災「1.17の記録」に阪神・淡路大震災の記録写真が公開されています。前回行ったようにこの写真にジオタグを追加して、その緯度経度にポイントを配置してみます。
上記サイトにおいて、画像データの一括ダウンロード＞CSVでダウンロードを行います。
こんな感じのCSVがダウンロードできました。緯度経度の項目があるにも関わらず、ここは空白になっています・・・仕方ないので以前行ったジオコーディングを使用したいと思います。ただ、場所が空白、もしくは、ジオコーディングできない（衣掛町5丁目2「こうべ地震災害対策広報」みたいな記述）写真は取込対象外とします。
主に以下のような技術を使用します。
①画像のダウンロード
②画像にジオタグ(geotag)の追加
③ジオコーディング
④ジオタグを追加した画像をポイントに変換
本エントリー作成時には「②画像にジオタグ(geotag)の追加」が正しく動作したのですが、現状、正しく動作が判明しないことがわかりました。以下エントリーの方法で同じことが実現できますので、「②画像にジオタグ(geotag)の追加」に関してはこちらを参照してください。
上述した「ポイントとなる技術」を盛り込んだサンプルコードです。パスなどはご自身の環境にあったものに書き換えてください。
# -*- coding: utf-8 -*- import pandas as pd import urllib2 import cookielib import os import pyexiv2 from PIL import Image import geocoder import arcpy def main(): #ダウンロードしたcsvの読込 df = pd.read_csv(r"D:\python\geotag\20160120142702.csv",encoding="SHIFT-JIS") #各種設定 folder = "D:\python\geotag" outFeatures = r"C:\ArcPySample\ArcPySample.gdb\geotag" badPhotosList = r"C:\ArcPySample\ArcPySample.gdb\empty" photoOption = "ONLY_GEOTAGGED" attachmentsOption = "" cj = cookielib.CookieJar() opener = urllib2.build_opener(urllib2.HTTPCookieProcessor(cj)) #CSVの行数分ループ for key,rowD in df.iterrows(): #リスエスト req = urllib2.Request(rowD.image) #出力先を指定 with open(os.path.join(folder,rowD.filename), 'wb') as f: f.write(opener.open(req).read()) #場所が無いデータは無視 if type(rowD.place) != float: #緯度経度取得 loc = get_coordinate(rowD.place) #井戸経度が取得できた場合 if loc.lat is not None: #ジオタグを画像に追加 set_gps_location(os.path.join("D:\python\geotag",rowD.filename),loc.lat,loc.lng) loc = None #ジオタグ付き写真をポイントに変換 arcpy.GeoTaggedPhotosToPoints_management(folder, outFeatures, badPhotosList, photoOption, attachmentsOption) def get_coordinate(location_name): try: #地名から座標を取得する ret = geocoder.google(location_name) except KeyError, e: print "KeyError(%s)" % str(e) return return ret def to_deg(value, loc): if value < 0: loc_value = loc[0] elif value > 0: loc_value = loc[1] else: loc_value = "" abs_value = abs(value) deg = int(abs_value) t1 = (abs_value-deg)*60 min = int(t1) sec = round((t1 - min)* 60, 5) return (deg, min, sec, loc_value) def set_gps_location(file_name, lat, lng): try: lat_deg = to_deg(lat, ["S", "N"]) lng_deg = to_deg(lng, ["W", "E"]) # 緯度、経度を10進法→60進法(度分秒)に変換 exiv_lat = (pyexiv2.Rational(lat_deg[0]*60+lat_deg[1],60),pyexiv2.Rational(lat_deg[2]*100,6000), pyexiv2.Rational(0, 1)) exiv_lng = (pyexiv2.Rational(lng_deg[0]*60+lng_deg[1],60),pyexiv2.Rational(lng_deg[2]*100,6000), pyexiv2.Rational(0, 1)) metadata = pyexiv2.ImageMetadata(file_name) metadata.read() #メタデータを付与 metadata["Exif.GPSInfo.GPSLatitude"] = exiv_lat metadata["Exif.GPSInfo.GPSLatitudeRef"] = lat_deg[3] metadata["Exif.GPSInfo.GPSLongitude"] = exiv_lng metadata["Exif.GPSInfo.GPSLongitudeRef"] = lng_deg[3] metadata["Exif.Image.GPSTag"] = 654 metadata["Exif.GPSInfo.GPSMapDatum"] = "WGS-84" metadata["Exif.GPSInfo.GPSVersionID"] = '2 0 0 0' #メタデータを上書き metadata.write() except: print u"エラー。たまに出る。とりあえず飛ばす" if __name__ == '__main__': main()
結果を確認してみます。それっぽい箇所にポイントが配置されましたが、全然関係ない場所にも配置されてしまっています。これは例えば「若松町」という場所が他の都道府県にも存在するという原因かと思います。完璧に行う場合はもう少し工夫が必要ですがとりあえず気にせず進みます（笑）
うまくいったっぽい箇所をGoogleMapと比較して確認してみます。「ポートターミナル周辺」でジオコーディング、GoogleMapで検索した結果です。
大体あってますかね。ちなみにこんな写真です。ジオタグもちゃんと追加されていますね。
ArcMapでポイントの個別属性から画像を参照することができます。
これだけで終わるのもあれなので、ここで一工夫してみようと思います。まずフィールドの追加を行います。
arcpy.AddField_management("geotag", "HTML_Path", "String")
こうなったらOKです。
その後、以下のようにフィールド演算を行います。
codeblock = """ def setHTMLPath(Path): return "<img src='" + Path + "'" + "Width='300'>" """ expression = "setHTMLPath(!Path!)" arcpy.CalculateField_management("geotag", "HTML_Path", expression,"PYTHON_9.3",codeblock)
こうなったらOKです。
ArcMapでHTMLポップアップ機能を使用します。
任意のポイントを選択するとポップアップが起動します。属性にもっと詳細な情報を入れたらより充実したマップになりますね。
今回わかったことは
①ダウンロードしたCSVの情報がかなり大雑把。そのため、ジオコーディングに一工夫が必要。
②ジオタグ追加→ポイント変換のコンボはかなり面白い！GPSの精度がさらに向上すれば、こういった案件も増えてきそうな予感。
③写真を見るとあらためて震災の怖さというものを感じた。
もっとオープンデータが充実すればさらに精度の高いマップを作成することができると思うので、今後そこを期待したいですね。
以上、「【ArcMapとArcpy】オープンデータを使って自分だけの地図を作ろう！ ～阪神・淡路大震災編～ 」でした。
本エントリー作成時にはサンプルコードが正しく動作したのですが、現状、正しく動作が判明しないことがわかりました。以下エントリーの方法で同じことが実現できますので、こちらを参照してください。
さて、本日は「Pythonで画像ファイルにジオタグ(geotag)を追加しよう！」です。GPSを搭載したスマートフォンなどは、撮影時にその場所の緯度経度情報を写真データに埋め込むことが可能で、この位置情報を「ジオタグ」と呼びます。このジオタグですが、注意が必要です。自宅で撮った写真をブログやSNSで公開する場合、ジオタグ付きの写真をそのままネットで公開してしまうと、撮影場所の緯度経度、つまりあなたの住所が丸わかりになってしまいます。悪用されたら怖いですね・・・
pyexiv2、PILが必要なので、インストールをお願いします。
神戸市のオープンデータを使用してみました。
1ファイルだけを対象としたサンプルコードです。
# -*- coding: utf-8 -*- import pyexiv2 from PIL import Image def to_deg(value, loc): if value < 0: loc_value = loc[0] elif value > 0: loc_value = loc[1] else: loc_value = "" abs_value = abs(value) deg = int(abs_value) t1 = (abs_value-deg)*60 min = int(t1) sec = round((t1 - min)* 60, 5) return (deg, min, sec, loc_value) def set_gps_location(file_name, lat, lng): lat_deg = to_deg(lat, ["S", "N"]) lng_deg = to_deg(lng, ["W", "E"]) # 緯度、経度を10進法→60進法(度分秒)に変換 exiv_lat = (pyexiv2.Rational(lat_deg[0]*60+lat_deg[1],60),pyexiv2.Rational(lat_deg[2]*100,6000), pyexiv2.Rational(0, 1)) exiv_lng = (pyexiv2.Rational(lng_deg[0]*60+lng_deg[1],60),pyexiv2.Rational(lng_deg[2]*100,6000), pyexiv2.Rational(0, 1)) metadata = pyexiv2.ImageMetadata(file_name) metadata.read() #メタデータを付与 metadata["Exif.GPSInfo.GPSLatitude"] = exiv_lat metadata["Exif.GPSInfo.GPSLatitudeRef"] = lat_deg[3] metadata["Exif.GPSInfo.GPSLongitude"] = exiv_lng metadata["Exif.GPSInfo.GPSLongitudeRef"] = lng_deg[3] metadata["Exif.Image.GPSTag"] = 654 metadata["Exif.GPSInfo.GPSMapDatum"] = "WGS-84" metadata["Exif.GPSInfo.GPSVersionID"] = '2 0 0 0' #メタデータを上書き metadata.write() #ファイル名、緯度、経度を指定 set_gps_location(r"D:\python\geotag\a020.jpg", 34.6874422,135.1758424)
結果を見ると、指定した緯度経度が設定されていますね。
ためしにArcMap上でジオタグ付き写真→ポイントに変換してみました。それっぽい場所に配置されましたね。
ポイントをクリックすると写真のポップアップを起動することもできます。これはサイズ調整が必要ですが・・・
これをうまく活用すれば震災マップなども作成することができますね。
以上、簡単ではありますが、「Pythonで画像ファイルにジオタグ(geotag)を追加しよう！」でした。
今回はfactbookを使用してみようと思います。factbookとは世界各国に関する情報を年鑑形式でまとめたアメリカ合衆国中央情報局 (CIA) の年次刊行物のことです。こんな便利なものが公開されているのですね。
以下サイトからダウンロードします。JSONファイルです。
factbook/factbook.json: World Factbook Coun... - GitHub
JSONなのでこんな感じですね。項目一覧が無いので、欲しい情報がどれなのか一つ一つ項目を確認しなければならないのがちょっとしんどいですね。
今回はその国の人口における65歳以上の割合、普通出生率、GDP($)を取得したいと思います。普通出生率とは日本で一般的な合計特殊出生率ではなく、人口1000人あたりにおける出生数らしいです。まずはJSONのデータを取得するサンプルです。
# -*- coding: utf-8 -*- import json from os.path import join,relpath from glob import glob import pandas as pd #ダウンロードしたフォルダのパスを指定 path = r"D:\python\factbook\factbook.json-master" dataList = [] #各大陸のフォルダをループ for fs in [relpath(x,path) for x in glob(join(path,'*'))]: #フォルダの場合 if fs.find('.') == -1: folder = join(path , fs) #フォルダ内のjsonファイルをループ for js in [relpath(x,folder) for x in glob(join(folder,'*'))]: if js != "xx.json": #jsonファイルを開く with open(join(folder , js), "r") as f: #ロード jsonData = json.load(f) #ここから階層をくだっていく for name in jsonData.keys(): groupDict = jsonData[name] if name == "Government": #この下のキーでループ for name2 in groupDict.keys(): groupDict2 = groupDict[name2] if name2 == "Country name": #この下のキーでループ for name3 in groupDict2.keys(): groupDict3 = groupDict2[name3] if name3 == "conventional short form": #この下のキーでループ for name4 in groupDict3.keys(): groupDict4 = groupDict3[name4] if name4 == "text": counryName = groupDict4 + "," elif name == "People and Society": #この下のキーでループ for name2 in groupDict.keys(): groupDict2 = groupDict[name2] if name2 == "Age structure": #この下のキーでループ for name3 in groupDict2.keys(): groupDict3 = groupDict2[name3] #65歳以上の割合 if name3 == "65 years and over": #この下のキーでループ for name4 in groupDict3.keys(): groupDict4 = groupDict3[name4] if name4 == "text": overSixtyFive = groupDict4[0:groupDict4.find("%")] + "," #普通出生率 elif name2 == "Birth rate": for name3 in groupDict2.keys(): groupDict3 = groupDict2[name3] if name3 == "text": birthRate = groupDict3[0:groupDict3.find(" births")] + "," elif name == "Economy": #この下のキーでループ for name2 in groupDict.keys(): groupDict2 = groupDict[name2] #GDP($) if name2 == "GDP (official exchange rate)": #この下のキーでループ for name3 in groupDict2.keys(): groupDict3 = groupDict2[name3] if name3 == "text": if groupDict3.find(' million') != -1: GDP = str(float(groupDict3[1:groupDict3.find(" million")]) * 1000000) elif groupDict3.find(' billion') != -1: GDP = str(float(groupDict3[1:groupDict3.find(" billion")]) * 1000000000) elif groupDict3.find(' trillion') != -1: GDP = str(float(groupDict3[1:groupDict3.find(" trillion")]) * 1000000000000) else: GDP = "0" dataList.append([counryName,overSixtyFive,birthRate,GDP]) #データフレーム作成 df = pd.DataFrame(dataList) df.columns = [u"counryName",u"overSixtyFive",u"birthRate","GDP"]
以前使用した世界地図を用意します。
まず、このレイヤにカラム追加します。
arcpy.AddField_management("TM_WORLD_BORDERS-0.3", "overSF", "Double") arcpy.AddField_management("TM_WORLD_BORDERS-0.3", "birthRate", "Double") arcpy.AddField_management("TM_WORLD_BORDERS-0.3", "GDP", "Double")
こんな感じになったら成功です。
次にJSONから取得した値をフィーチャの属性に反映させます。
for key,rowD in df.iterrows(): cursorJ = arcpy.UpdateCursor("TM_WORLD_BORDERS-0.3") for rowJ in cursorJ: #国名が微妙に異なるのはとりあえず無視(GambiaとThe Gambiaとか) if rowD.counryName == rowJ.NAME: rowJ.setValue("overSF", rowD.overSixtyFive) rowJ.setValue("birthRate", rowD.birthRate) rowJ.setValue("GDP", rowD.GDP) cursorJ.updateRow(rowJ)
こんな感じになったら成功です。
データもそろったので分析をしてみたいと思います。
レイヤの設定を以下のようにします。
やはりわが日本が厳しい状態ですね。ヨーロッパ、カナダなどの先進国は高齢化していて、その反面、アフリカや南アジアなどの新興国は高齢者の割合が低いことが分かりました（アフリカなどは平均寿命の関係もありそうですが）。
前述した通り、普通出生率とは日本で一般的な合計特殊出生率ではなく、人口1000人あたりにおける出生数らしいです。
レイヤの設定を以下のようにします。
これもわが日本を含め先進国は厳しい状態ですね。
レイヤの設定を以下のようにします。
GDPはアメリカと中国が突出していて、次に日本を含めた先進国が高いという結果が出ました。まぁこれはみなさんよくご存知かもしれませんね。また、アフリカ、東欧、東南アジアに今後の伸びしろがあるのがわかります。東南アジア経済は今熱いですしね！
今回は65歳以上の割合、普通出生率、GDP($)のデータを使用しましたが、factbookにはその他にもたくさんのデータがあります。例えば、エネルギー、軍事、経済、政治などですね。これらをうまく活用すれば世界を色々な視点から分析することが可能になりますね。
以上、「【ArcMapとArcpy】オープンデータを使って自分だけの地図を作ろう！ ～factbook編～」でした。
本日はPythonでのArcobjectsの使用方法を書いてみたいと思います。PythonだとArcpyしか使えず、ちょっと難しいことをやろうと思うと途端にその貧弱さが露呈してしまいましたが、ArcObjectsを使えば大丈夫？！なはずです。それじゃあ最初から.NETで開発しろよというツッコミはなしで（我らがPythonを使いたいのです！）。
ArcGISのバージョンごとにこんなものが用意されています。
・9.3, 9.4 - http://pierssen.com/arcgis/upload/misc/python_arcobjects.zip
・10.0 - http://www.pierssen.com/arcgis10/upload/python/snippets100.py
・10.1 - http://www.pierssen.com/arcgis10/upload/python/snippets101.py
・10.2 - http://www.pierssen.com/arcgis10/upload/python/snippets102.py
これらのうちいずれかを例えば、C:\Python27\ArcGIS10.2\Scriptsに配置してください。# -*- coding: utf-8 -*- を忘れずに！
https://pypi.python.org/pypi/comtypesにソースがありますので、そこからpython setup.py install からのeasy_install comtypes-1.1.2-py2.7.egg ですね。
ポイントを作ってその座標を設定するサンプルです。
# -*- coding: utf-8 -*- from comtypes.client import GetModule, CreateObject from snippets102 import GetLibPath, InitStandalone m = GetModule(GetLibPath() + "esriGeometry.olb") p = CreateObject(m.Point, interface=m.IPoint) p.PutCoords(2,3) print p.X, p.Y
詳細は使用方法についてはこれから調査してみたいと思います。
以上、本日は「PythonでArcObjectsを使ってみよう！」でした。
今回はオープンデータについてです。オープンデータとは誰もが自由に利用でき、再利用や再配布が許可されているデータのことを指します。つまり、「著作権フリーで自由に利用可能」ということです。素晴らしいですね！オープンデータをうまく活用すれば地図上にこんなことをすることが可能です。
①「世界のサッカー場を配置してみる」
②「世界のダイビングポイントを配置してみる」
③「世界のお酒の名店を配置してみる」
ただし、それぞれ座標を調べる必要があるので、そういったものが公開されているデータを利用するのが現実的かと思います。←①を作りたいのですが、公開されているデータがなく断念しました。
とりあえず、今回は世界の空港が配置された世界地図を作ってみようと思います。まずは素の世界地図をthematicmapping.orgからダウンロードします。
いい感じですね！
続いて、空港のデータをOpenFlights: Airport and airline dataからダウンロードします。
datになっていますね。
コピペしてCSVを作成します。以下のような感じですね。
そして、マップレイヤに作成したCSVを追加してください。
CSVを右クリックして「XYデータの表示」をクリックしてください。←これArcpyでできないのかなぁ。
こんな感じでOKをおしてください。
エラーが出ますが、無視してかまいません。
ポイントが配置されました。空港ってこんなあるんですね（笑）
日本だけでもこんなにあります。
とりあえず国別の空港数と1空港あたりの利用人数を抽出してみました。
numpyとpandasが必要です。
# -*- coding: utf-8 -*- import arcpy import pandas as pd import numpy as np #フィーチャの属性をpandasに食べさせる arrAir = arcpy.da.TableToNumPyArray("AirPort", ('Country')) df = pd.DataFrame(arrAir) #国ごとの空港数を集計する dfAir = pd.DataFrame(df['Country'].value_counts()) dfAir.columns = ['AirPortNum'] dfAir['Country'] = 0 for key,rowD in dfAir.iterrows(): dfAir['Country'][key] = key #フィーチャの属性をpandasに食べさせる arrCountry = arcpy.da.TableToNumPyArray("TM_WORLD_BORDERS-0.3", ('NAME','POP2005')) dfCountry = pd.DataFrame(arrCountry) dfCountry.columns = ['Country', 'POP2005'] #データフレームを結合させる dfMerge= pd.merge(dfAir,dfCountry,on='Country',how='right') #1空港あたりの利用人数を計算 dfMerge['CAL'] = dfMerge['POP2005'] / dfMerge['AirPortNum'] #CSV出力 dfMerge.to_csv('D:\python\global\AirPortByCountry.csv',encoding="SHIFT-JIS")
出力したCSVを確認してみるとアメリカが圧倒的ですね。AirPortNumが空港数、POP2005が2005年の人口、CALが1空港あたりの利用人数です。これを元に国の色分けをしたらおもしろそうなので、以下スクリプトを追加してください。
#カラム追加 arcpy.AddField_management("TM_WORLD_BORDERS-0.3", "CAL", "Double") for key,rowD in dfMerge.iterrows(): cursorJ = arcpy.UpdateCursor("TM_WORLD_BORDERS-0.3") for rowJ in cursorJ: if rowD.Country == rowJ.NAME: rowJ.setValue("CAL", rowD.CAL) cursorJ.updateRow(rowJ)
こうなったら成功です。
シンボル設定をこのようにしてみました。
結果を見ると中国、インドなど人口の多い国が1空港あたりの利用人数が多いということがわかりました。また、欧米は1空港あたりの利用人数が少ないことがわかりました。なんでもヨーロッパなんかは線路を引くのが面倒だから飛行機を飛ばしちまえみたいな感覚だとかなんとか。
データ＞データのエクスポート を行います。これでデータの作成が完了しました。
こんな感じでシンボルを設定してもおもしろいですね。
上記ではShapeに出力しましたが、GDBに保存して、シンボル設定をしたレイヤファイルを作成した方がいいかもしれませんね。そして、今回はたまたま空港のデータを使用しましたが、色々探せばおもしろいデータを発見することができそうですね。例えば、国ごとの犯罪率や地震発生率などを視覚化することができますし、また、サッカーW杯優勝回数やその国に属するクラブの世界ランキングを考慮した上の国ごとの世界ランキングなんかも作れると思います。その場合、やはりヨーロッパ、南米が強いということが視覚的にわかるのではないでしょうか。次回以降も色々やってみたいと思います。
以上、本日は「【ArcMapとArcpy】オープンデータを使って自分だけの地図を作ろう！ ～世界の空港編～」でした。
