読者です 読者をやめる 読者になる 読者になる

GIS奮闘記

現役GISエンジニアの技術紹介ブログ。主にPythonを使用。

Arcpyを使用してJ-SHISの表層地盤データを分析してみよう(横須賀編) ~ArcpyでShape読込、Shape出力、属性検索、空間検索、フィールド演算~

本日はJ-SHISと地盤についてです。ご存じの方も多いかも知れませんが、J-SHISとは「地震防災に資することを目的に、日本全国の「地震ハザードの共通情報 基盤」として活用されることを目指して作られたサービス」です。ここから全国の地盤データをダウンロードすることができます。詳細は以下サイトをご参照ください。
http://www.j-shis.bosai.go.jp/map/JSHIS2/download.html?lang=jp


横須賀を含むこの区域の表層地盤データをダウンロードします。
f:id:sanvarie:20160106110052p:plain


ただし、これだけだと横須賀以外の地域が含まれてしまいます。ですので、横須賀の表層地盤データのみを抽出したいと思います。

まず、ESRIから全国地図をダウンロードします。
全国市区町村界データ | 製品 | ESRIジャパン


以下pythonスクリプトをArcMap上で実行してください。ダウンロードした表層地盤データと全国地図を使用して横須賀の表層地盤データを抽出します。Arcpyを用いてこのようにけっこう色々なことをします。
・Shapeファイルの読み込み
・属性検索
・空間検索
・Shape出力
・読込したShapeファイルの削除
・フィールド演算

# -*- coding: utf-8 -*-
import arcpy

#J-SHISのデータ(表層地盤)
inLyr1 = ur"D:\python\地盤高\Z-V3-JAPAN-AMP-VS400_M250-SHAPE-5239.shp"

#ESRIのデータ(日本地図)
inLyr2 = ur"D:\python\地盤高\japan_ver80.shp"

mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd)[0]
addLayer1 = arcpy.mapping.Layer(inLyr1)
addLayer2 = arcpy.mapping.Layer(inLyr2)

#Shapeファイルの読込
arcpy.mapping.AddLayer(df, addLayer1, "BOTTOM")
arcpy.mapping.AddLayer(df, addLayer2, "BOTTOM")

#属性検索(横須賀市を検索)
zokuseiSearch = arcpy.SelectLayerByAttribute_management(addLayer2.name,"NEW_SELECTION","SIKUCHOSON = '横須賀市'")

#空間検索(横須賀に交差する地盤高データを検索)
kukanSearch = arcpy.SelectLayerByLocation_management(addLayer1.name,"INTERSECT",zokuseiSearch,"","NEW_SELECTION")

#空間検索した結果をShape出力
arcpy.FeatureClassToFeatureClass_conversion(kukanSearch,
                                            u"D:\python\地盤高",
                                            "zibandaka")

#出力したShape以外のレイヤを削除
for df in arcpy.mapping.ListDataFrames(mxd):
    for lyr in arcpy.mapping.ListLayers(mxd, "", df):
        if lyr.name != "zibandaka":
            arcpy.mapping.RemoveLayer(df, lyr)

#コードブロック
codeblock = """
def GetBichikei(JCODE):
    if JCODE == 1:
        return "山地"
    elif JCODE == 2:
        return "山麓地"
    elif JCODE == 3:
        return "丘陵"
    elif JCODE == 4:
        return "火山地"
    elif JCODE == 5:
        return "火山山麓地"
    elif JCODE == 6:
        return "火山性丘陵"
    elif JCODE == 7:
        return "岩石台地"
    elif JCODE == 8:
        return "砂礫質台地"
    elif JCODE == 9:
        return "ローム台地"
    elif JCODE == 10:
        return "谷底低地"
    elif JCODE == 11:
        return "扇状地"
    elif JCODE == 12:
        return "自然堤防"
    elif JCODE == 13:
        return "後背湿地"
    elif JCODE == 14:
        return "旧河道"
    elif JCODE == 15:
        return "三角州・海岸低地"
    elif JCODE == 16:
        return "砂州・砂礫州"
    elif JCODE == 17:
        return "砂丘"
    elif JCODE == 18:
        return "砂州・砂丘間低地"
    elif JCODE == 19:
        return "干拓地"
    elif JCODE == 20:
        return "埋立地"
    elif JCODE == 21:
        return "磯・岩礁"
    elif JCODE == 22:
        return "河原"
    elif JCODE == 23:
        return "河道"
    elif JCODE == 24:
        return "湖沼"
    else:
        return "不明"
"""

# 条件式を設定
expression = "GetBichikei(!JCODE!)"

#微地形分類名カラムを追加
arcpy.AddField_management("zibandaka", "bichikei", "TEXT")

#フィールド演算で微地形分類名を設定
arcpy.CalculateField_management("zibandaka", "bichikei", expression, "PYTHON_9.3", codeblock)

以下のような結果が得られました。

f:id:sanvarie:20160106135022p:plain

日本地図と重ねあわせるとこんな感じです。ぴったり横須賀の表層地盤のみ出力できていますね。

f:id:sanvarie:20160106131731p:plain

上記プログラムで微地形分類コードの名称を取得していますが、微地形分類とは簡単にいうと、地盤の種類ですね。詳細は以下サイトで確認できます。
http://www.j-shis.bosai.go.jp/map/JSHIS2/data/DOC/DataFileRule/A-RULES.pdf

次は以下のようにArcMapでzibandakaレイヤのシンボルのプロパティを変更してください。どうやらレイヤのsymbologyTypeがOTHERSになっていると、Arcpyではレンダリングの編集はできないみたいです。なんて貧弱なんだ!←ESRIさんお願いしますよ。

f:id:sanvarie:20160106135322p:plain

このように横須賀の表層地盤の視覚化に成功しました(なぜ横須賀を選んだのかは秘密です)!これを見ると横須賀は丘陵、沿岸部に埋立地が多いことがわかりました。ちなみにこの微地形分類は震災時の被害予測データとして使われたりします。また、これに道路や建物のデータを重ねあわせると視覚的にはさらに面白くなりそうですね。

f:id:sanvarie:20160106154433p:plain

以上、「Arcpyを使用してJ-SHISの表層地盤データを分析してみよう(横須賀編) ~ArcpyでShape読込、Shape出力、属性検索、空間検索、フィールド演算~」でした。