砂防基盤図の古い航空写真を刷新

砂防基盤図の古い航空写真を刷新：QGISで地理院シームレス写真を1/1,250図郭ごとに0.25 m/pxで一括生成し、土石流区域設定システムへ取り込む

砂防基盤図の古い航空写真を刷新：QGISで地理院シームレス写真を1/1,250図郭ごとに0.25 m/pxで一括生成し、土石流区域設定システムへ取り込む

目的

既存の砂防基盤図に付属する航空写真が古いため、国土地理院「シームレス写真」を用いて最新の画像を国土基本図郭（1/1,250）ごとに作成し、砂防フロンティア「土石流区域設定システム」へ取り込むための素材（JPG＋JGW または GeoTIFF＋TFW）を一括出力します。

成果物

  • 各図郭ごとに ort_MAPID.jpg（または .tif）と JGW/TFW を生成
  • 解像度 0.25 m/px 固定：1/1,250 図郭（約 1000×750 m）は4000×3000 pxで統一
  • 出力先は図郭シェープ（.shp）と同じフォルダ

前提

  • QGIS 3.x
  • 図郭ポリゴン（例：ZUKAKU_1250_7）CRS はメートル系（平面直角）
  • 地理院シームレス写真（XYZ）URL：https://cyberjapandata.gsi.go.jp/xyz/seamlessphoto/{z}/{x}/{y}.jpg

手順（概要）

  1. 図郭シェープをQGISへ追加（必要なら対象図郭のみ選択）。
  2. 下記の PyQGIS スクリプトを「Python コンソール → エディタ」に貼り付けて F5 実行。
  3. 図郭 .shp と同じフォルダに ort_*.jpg＋.jgw（または .tif＋.tfw）が出力されます。

PyQGIS スクリプト（0.25 m/px 固定）

# 1/1,250 図郭に合わせて 0.25 m/px で切り出し（JPG + JGW）
# 1/2,500 図郭も同じ m/px で書き出し可能（8000×6000 px）

import os
from qgis.core import (
    QgsProject, QgsRasterLayer, QgsWkbTypes, QgsUnitTypes,
    QgsMapSettings, QgsMapRendererParallelJob
)
from qgis.PyQt.QtCore import QSize

# ===== 設定 =====
POLY_NAME   = 'ZUKAKU_1250_7'                 # 図郭レイヤ名
NAME_FIELD  = 'MAPID'                          # 出力名に使う属性（例：HD4431）
PIXEL_M     = 0.25                             # ★ 0.25 m/px 固定
XYZ_NAME    = 'GSI_seamless'                   # 既に追加済みならこの名前で取得
XYZ_URI     = ('type=xyz&'
               'url=https://cyberjapandata.gsi.go.jp/xyz/seamlessphoto/{z}/{x}/{y}.jpg&'
               'zmax=18&zmin=14')
OUTPUT_FMT  = 'JPG'                            # 'JPG' または 'TIFF'
JPEG_QUALITY= 95                               # JPG の画質
SET_DPI     = 254                              # 印刷用のDPIメタ（任意・見え方は不変）
# ===============

poly = QgsProject.instance().mapLayersByName(POLY_NAME)[0]
if poly.crs().mapUnits() != QgsUnitTypes.DistanceMeters:
    raise ValueError('図郭レイヤのCRSはメートル系（平面直角など）にしてください')

out_dir = os.path.dirname(poly.source().split('|')[0])
print(f'出力フォルダ: {out_dir}')

lst = QgsProject.instance().mapLayersByName(XYZ_NAME)
xyz = lst[0] if lst else QgsRasterLayer(XYZ_URI, XYZ_NAME, 'wms')
if not lst:
    if not xyz.isValid(): raise ValueError('シームレス写真レイヤ作成に失敗')
    QgsProject.instance().addMapLayer(xyz, False)

feats = poly.selectedFeatures() or list(poly.getFeatures())
print(f'対象: {len(feats)} 件  |  解像度: {PIXEL_M:.2f} m/px')

for f in feats:
    code = str(f[NAME_FIELD] or f.id()).strip()
    ext  = f.geometry().boundingBox()
    w = int(round(ext.width()  / PIXEL_M))
    h = int(round(ext.height() / PIXEL_M))

    ms = QgsMapSettings()
    ms.setLayers([xyz])
    ms.setDestinationCrs(poly.crs())
    ms.setExtent(ext)
    ms.setOutputSize(QSize(w, h))

    job = QgsMapRendererParallelJob(ms)
    job.start(); job.waitForFinished()
    img = job.renderedImage()

    if SET_DPI:
        dpm = int(SET_DPI / 0.0254)
        img.setDotsPerMeterX(dpm); img.setDotsPerMeterY(dpm)

    stem = f'ort_{code}'
    if OUTPUT_FMT.upper() == 'TIFF':
        out_path = os.path.join(out_dir, stem + '.tif')
        ok = img.save(out_path, 'TIFF')
        wld = os.path.join(out_dir, stem + '.tfw')
    else:
        out_path = os.path.join(out_dir, stem + '.jpg')
        ok = img.save(out_path, 'JPG', JPEG_QUALITY)
        wld = os.path.join(out_dir, stem + '.jgw')

    xres = PIXEL_M; yres = -PIXEL_M
    xul  = ext.xMinimum() + xres/2.0
    yul  = ext.yMaximum() - PIXEL_M/2.0
    with open(wld, 'w', encoding='ascii') as fp:
        fp.write(f'{xres:.10f}\n0.0\n0.0\n{yres:.10f}\n{xul:.10f}\n{yul:.10f}\n')

    print(f'✓ {code} → {os.path.basename(out_path)} [{w}×{h}px, {PIXEL_M:.2f} m/px]')

print('完了')
取り込み（砂防フロンティア・土石流区域設定システム）に向けたポイント
・システム側が受け付けるラスタは多くの場合「JPG+JGW」または「GeoTIFF(+TFW)」。本記事の出力は両対応です。
・CRS は作業系（例：平面直角 第○系）で統一してください。システム側要求と合わせるのが無難です。
・図郭ごとにファイル名へ MAPID を付けておくと、位置の突合や登録時の検索が容易になります。
・既存の砂防基盤図に重ねて目視確認し、欠落や雲の多いタイルがあれば該当図郭だけ再出力します。

運用ヒント

  • より精細に：PIXEL_M=0.20（PC負荷に注意）。逆に軽く：0.30–0.40
  • 一部図郭のみ更新：QGISで対象図郭を選択して実行すると、選択分だけ出力
  • 画像メタの DPI は印刷用。地図の精細さは JGW/TFW の 1行目（m/px）で決まります。

ライセンス表記

出典：地理院タイル（シームレス写真）。利用規約に従い、出力物・画面等に出典を明記してください。

コメント

  1. eon2025年9月9日 19:00

    下記は、出力TIFバージョン


    # 1/1,250 図郭に合わせて 0.25 m/px で切り出し（TIFF + TFW）
    # 1/2,500 図郭も同じ m/px で書き出し可能（8000×6000 px）

    import os
    from qgis.core import (
    QgsProject, QgsRasterLayer, QgsWkbTypes, QgsUnitTypes,
    QgsMapSettings, QgsMapRendererParallelJob
    )
    from qgis.PyQt.QtCore import QSize

    # ===== 設定 =====
    POLY_NAME = 'ZUKAKU_1250_7' # 図郭レイヤ名
    NAME_FIELD = 'MAPID' # 出力名に使う属性（例：HD4431）
    PIXEL_M = 0.25 # ★ 0.25 m/px 固定
    XYZ_NAME = 'GSI_seamless' # 既に追加済みならこの名前で取得
    XYZ_URI = ('type=xyz&'
    'url=https://cyberjapandata.gsi.go.jp/xyz/seamlessphoto/{z}/{x}/{y}.jpg&'
    'zmax=18&zmin=14')
    SET_DPI = 254 # 印刷用のDPIメタ（任意・見え方は不変）
    # ===============

    poly = QgsProject.instance().mapLayersByName(POLY_NAME)[0]
    if poly.crs().mapUnits() != QgsUnitTypes.DistanceMeters:
    raise ValueError('図郭レイヤのCRSはメートル系（平面直角など）にしてください')

    out_dir = os.path.dirname(poly.source().split('|')[0])
    print(f'出力フォルダ: {out_dir}')

    lst = QgsProject.instance().mapLayersByName(XYZ_NAME)
    xyz = lst[0] if lst else QgsRasterLayer(XYZ_URI, XYZ_NAME, 'wms')
    if not lst:
    if not xyz.isValid(): raise ValueError('シームレス写真レイヤ作成に失敗')
    QgsProject.instance().addMapLayer(xyz, False)

    feats = poly.selectedFeatures() or list(poly.getFeatures())
    print(f'対象: {len(feats)} 件 | 解像度: {PIXEL_M:.2f} m/px')

    for f in feats:
    code = str(f[NAME_FIELD] or f.id()).strip()
    ext = f.geometry().boundingBox()
    w = int(round(ext.width() / PIXEL_M))
    h = int(round(ext.height() / PIXEL_M))

    ms = QgsMapSettings()
    ms.setLayers([xyz])
    ms.setDestinationCrs(poly.crs())
    ms.setExtent(ext)
    ms.setOutputSize(QSize(w, h))

    job = QgsMapRendererParallelJob(ms)
    job.start(); job.waitForFinished()
    img = job.renderedImage()

    if SET_DPI:
    dpm = int(SET_DPI / 0.0254)
    img.setDotsPerMeterX(dpm); img.setDotsPerMeterY(dpm)

    stem = f'ort_{code}'
    out_path = os.path.join(out_dir, stem + '.tif')
    ok = img.save(out_path, 'TIFF')
    if not ok:
    raise IOError(f'TIFF保存に失敗: {out_path}')

    # ワールドファイル（TFW）
    wld = os.path.join(out_dir, stem + '.tfw')
    xres = PIXEL_M; yres = -PIXEL_M
    xul = ext.xMinimum() + xres/2.0
    yul = ext.yMaximum() - PIXEL_M/2.0
    with open(wld, 'w', encoding='ascii') as fp:
    fp.write(f'{xres:.10f}\n0.0\n0.0\n{yres:.10f}\n{xul:.10f}\n{yul:.10f}\n')

    print(f'✓ {code} → {os.path.basename(out_path)} [{w}×{h}px, {PIXEL_M:.2f} m/px]')

    print('完了')

    返信削除
  2. eon2025年11月10日 23:30

    QGISで、出力したい図郭レイヤ（toyama_lvl2500_epsg6675 など）をアクティブにする（レイヤパネルでクリックして青選択）。

    そのレイヤ上で必要な図郭だけ選択。

    下記コードを Python コンソールで実行。

    返信削除
  3. eon2025年11月10日 23:31

    import os
    from qgis.core import (
    QgsProject,
    QgsRasterLayer,
    QgsUnitTypes,
    QgsMapSettings,
    QgsMapRendererParallelJob,
    )
    from qgis.PyQt.QtCore import QSize

    # ===== 設定 =====
    NAME_FIELD = 'code2500' # ファイル名に使う図郭コード列
    PIXEL_M = 0.25
    XYZ_NAME = 'GSI_seamless'
    XYZ_URI = (
    'type=xyz&'
    'url=https://cyberjapandata.gsi.go.jp/xyz/seamlessphoto/{z}/{x}/{y}.jpg&'
    'zmax=18&zmin=14'
    )
    SET_DPI = 254
    # ===============

    # アクティブな図郭レイヤを取得
    poly = iface.activeLayer()
    if not poly:
    raise ValueError('図郭レイヤをアクティブにしてから実行してください')

    print(f'使用レイヤ: {poly.name()}')

    # 座標単位チェック
    if poly.crs().mapUnits() != QgsUnitTypes.DistanceMeters:
    raise ValueError('図郭レイヤのCRSは平面直角などメートル系である必要があります')

    # 出力フォルダ = 図郭レイヤのソースと同じ場所
    src_path = poly.source().split('|')[0]
    out_dir = os.path.dirname(src_path)
    if not out_dir:
    out_dir = os.path.expanduser('~')
    print(f'出力フォルダ: {out_dir}')

    # GSIシームレスレイヤ取得
    xyz_layers = QgsProject.instance().mapLayersByName(XYZ_NAME)
    if xyz_layers:
    xyz = xyz_layers[0]
    else:
    xyz = QgsRasterLayer(XYZ_URI, XYZ_NAME, 'wms')
    if not xyz.isValid():
    raise ValueError('GSIシームレス写真レイヤの作成に失敗しました')
    QgsProject.instance().addMapLayer(xyz, False)

    # 対象図郭 = 選択されていれば選択のみ / なければ全件
    selected = poly.selectedFeatures()
    if selected:
    feats = selected
    print(f'選択図郭のみ出力します: {len(feats)} 件')
    else:
    feats = list(poly.getFeatures())
    print(f'※選択なしのため全図郭 {len(feats)} 件を対象とします')

    print(f'解像度: {PIXEL_M:.2f} m/px')

    for f in feats:
    code = str(f[NAME_FIELD] or f.id()).strip()
    if not code:
    continue

    ext = f.geometry().boundingBox()
    w = int(round(ext.width() / PIXEL_M))
    h = int(round(ext.height() / PIXEL_M))

    ms = QgsMapSettings()
    ms.setLayers([xyz])
    ms.setDestinationCrs(poly.crs())
    ms.setExtent(ext)
    ms.setOutputSize(QSize(w, h))

    job = QgsMapRendererParallelJob(ms)
    job.start()
    job.waitForFinished()
    img = job.renderedImage()

    if SET_DPI:
    dpm = int(SET_DPI / 0.0254)
    img.setDotsPerMeterX(dpm)
    img.setDotsPerMeterY(dpm)

    stem = f'ort_{code}'
    out_path = os.path.join(out_dir, stem + '.tif')
    if not img.save(out_path, 'TIFF'):
    raise IOError(f'TIFF保存に失敗: {out_path}')

    # ワールドファイル（TFW）
    wld = os.path.join(out_dir, stem + '.tfw')
    xres = PIXEL_M
    yres = -PIXEL_M
    xul = ext.xMinimum() + xres / 2.0
    yul = ext.yMaximum() - PIXEL_M / 2.0
    with open(wld, 'w', encoding='ascii') as fp:
    fp.write(f'{xres:.10f}\n')
    fp.write('0.0\n')
    fp.write('0.0\n')
    fp.write(f'{yres:.10f}\n')
    fp.write(f'{xul:.10f}\n')
    fp.write(f'{yul:.10f}\n')

    print(f'✓ {code} → {os.path.basename(out_path)} [{w}×{h}px]')

    print('完了')

    返信削除
  4. eon2025年11月11日 19:26

    import os
    from qgis.core import (
    QgsProject, QgsRasterLayer, QgsUnitTypes,
    QgsMapSettings, QgsMapRendererParallelJob
    )
    from qgis.PyQt.QtCore import QSize

    # ===== 設定 =====
    NAME_FIELD = 'code2500' # 図郭コード列名（例: GE881 等）
    PIXEL_M = 0.25 # m/px
    XYZ_NAME = 'GSI_seamless'
    XYZ_URI = (
    'type=xyz&'
    'url=https://cyberjapandata.gsi.go.jp/xyz/seamlessphoto/{z}/{x}/{y}.jpg&'
    'zmax=18&zmin=14&'
    'crs=EPSG:3857' # ★ 明示
    )
    SET_DPI = 254
    # ===============

    # --- 図郭レイヤ取得 ---
    poly = iface.activeLayer()
    if not poly:
    raise ValueError('図郭レイヤをアクティブにしてから実行してください')
    print(f'使用レイヤ: {poly.name()}')

    # --- CRS 単位チェック（m 必須）---
    if poly.crs().mapUnits() != QgsUnitTypes.DistanceMeters:
    raise ValueError('図郭レイヤのCRSはメートル系（平面直角など）が必要です')

    # --- 出力フォルダ決定 ---
    src_path = poly.source().split('|')[0]
    out_dir = os.path.dirname(src_path) if src_path else ''
    if not out_dir:
    out_dir = os.path.expanduser('~')
    if not os.path.isdir(out_dir) or not os.access(out_dir, os.W_OK):
    raise IOError(f'出力フォルダに書き込みできません: {out_dir}')
    print(f'出力フォルダ: {out_dir}')

    # --- XYZ（GSIシームレス） ---
    xyz_layers = QgsProject.instance().mapLayersByName(XYZ_NAME)
    if xyz_layers:
    xyz = xyz_layers[0]
    else:
    xyz = QgsRasterLayer(XYZ_URI, XYZ_NAME, 'wms') # provider=wms で OK（type=xyz）
    if not xyz.isValid():
    raise ValueError('GSIシームレス写真レイヤの作成に失敗しました（ネットワーク/URIを確認）')
    QgsProject.instance().addMapLayer(xyz, False)

    # --- 対象図郭一覧 ---
    selected = poly.selectedFeatures()
    feats = selected if selected else list(poly.getFeatures())
    print(('選択図郭のみ' if selected else '全図郭') + f'を出力します: {len(feats)} 件')
    print(f'解像度: {PIXEL_M:.2f} m/px')

    # --- ループ処理 ---
    for f in feats:
    # 図郭コード
    try:
    code_val = f[NAME_FIELD]
    except Exception:
    raise KeyError(f'属性「{NAME_FIELD}」が見つかりません。列名を確認してください。')

    code = str(code_val if code_val not in (None, '') else f.id()).strip()
    if not code:
    print(' → 属性が空のためスキップ')
    continue

    # 図郭範囲
    ext = f.geometry().boundingBox()
    # まれに幅/高さ 0 になるのを防ぐ
    if ext.width() == 0 or ext.height() == 0:
    ext.grow(0.01) # 1cm バッファ
    w = int(round(ext.width() / PIXEL_M))
    h = int(round(ext.height() / PIXEL_M))
    if w <= 0 or h <= 0:
    raise ValueError(f'ピクセルサイズが0または負です (code={code}, w={w}, h={h}, ext={ext.toString()})')

    # レンダリング設定
    ms = QgsMapSettings()
    ms.setLayers([xyz])
    ms.setDestinationCrs(poly.crs()) # 図郭CRSへ再投影して出力
    ms.setExtent(ext)
    ms.setOutputSize(QSize(w, h))

    # レンダリング
    job = QgsMapRendererParallelJob(ms)
    job.start(); job.waitForFinished()
    img = job.renderedImage()

    # DPI（任意）
    if SET_DPI:
    dpm = int(SET_DPI / 0.0254)
    img.setDotsPerMeterX(dpm)
    img.setDotsPerMeterY(dpm)

    # 保存
    stem = f'ort_{code}'
    out_path = os.path.join(out_dir, stem + '.tif')
    ok = img.save(out_path, 'TIFF')
    if not ok:
    raise IOError(f'TIFF保存に失敗: {out_path}')

    # ワールドファイル（TFW）
    wld = os.path.join(out_dir, stem + '.tfw')
    xres = PIXEL_M
    yres = -PIXEL_M
    xul = ext.xMinimum() + xres / 2.0
    yul = ext.yMaximum() - PIXEL_M / 2.0
    with open(wld, 'w', encoding='ascii') as fp:
    fp.write(f'{xres:.10f}\n0.0\n0.0\n{yres:.10f}\n{xul:.10f}\n{yul:.10f}\n')

    print(f'✓ {code} → {os.path.basename(out_path)} [{w}×{h}px]')

    print('完了')

    返信削除

コメントを投稿

このブログの人気の投稿

石川県：土砂災害（特別）警戒区域＋ＣＳ立体図

石川県の土砂災害（特別）警戒区域まとめ 2025 年 5 月 16 日 時点 石川県砂防課が公表する 土砂災害（特別）警戒区域指定状況 （PDF）を整理し、市町別に 土石流・地すべり・急傾斜地崩壊 の警戒区域（イエロー）と 特別警戒区域（レッド）を集計しました。 市町 警戒区域 (イエロー) 特別警戒区域 (レッド) 合計 土石流 地すべり 急傾斜 土石流 地すべり 急傾斜 輪島市 374 217 375 313 0 372 966 金沢市 230 82 247 170 0 235 559 能登町 240 68 218 210 0 209 526 七尾市 192 45 201 165 0 201 438 珠洲市 91 76 164 83 0 163 331 小松市 181 13 66 162 0 64 260 白山市 139 20 91 73 0 84 250 志賀町 130 15 100 106 0 99 245 穴水町 93 12 115 79 0 113 220 加賀市 78 21 101 65 0 96 200 津幡町 78 48 57 56 0 56 183 中能登町 130 4 36 99 0 35 170 能美市 64 1 41 59 0 40 106 宝達志水町 33 30 53 26 0 52 116 かほく市 12 2 30 5 0 30 44 川北町 0 0 0 0 0 0 0 野々市市 0 0 0 0 0 0 0 内灘町 0 0 5 0 0 5 5 合計 2065 654 1900 1671 0 1854 4619 ＣＳ立体図＋土砂災害警戒区域 MAP を開く 関連資料・データ ▸ 土砂災害（特別）警戒区域指定状況 PDF（R7.5.16） ▸ 石川県土砂災害（特別）警戒区域図 オープンデータ（SHP・GeoJSON 等） ▸ CS 立体図を使った 地形判読マニュアル（林野庁） ハザードマップの入手方法 お住まいの 市役所・町役場の防災担当窓口 や自...
続きを読む

日本でよく使う EPSG コード 一覧

1. 経緯度（地理）座標系 区分 名称 EPSG 主な用途・備考 世界標準（地理） WGS 84（緯度・経度） 4326 GNSS/GPS の測位やデータ交換で使用。単位は「度」。 Web地図（投影） WGS 84 / Pseudo-Mercator（Web Mercator） 3857 Google/OSM/国土地理院タイルなどの 表示 に使用。単位は「メートル」。別名：Popular Visualization / Google Mercator。 日本最新 JGD2011（緯度・経度） 6668 測地成果2011。公共測量の基準。平面直角座標系は下表参照。 過去の全国基準 JGD2000（緯度・経度） 4612 2002–2011 年に使用。JGD2011 と混在時は変換推奨。 旧日本測地系 Tokyo Datum 4301 1940年代〜2002年まで使用。WGS84とは数百mのズレあり。 2. 平面直角座標系（Gauss–Krüger 投影） 系 (ローマ数字) 系 (数字) 中央経線 JGD2011 EPSG JGD2000 EPSG 適用都道府県（概要） I 1 129°30′E 6669 2443 長崎、鹿児島西部 II 2 131°00′E 6670 2444 福岡・佐賀・熊本・大分・宮崎、鹿児島...
続きを読む

設計定数を求めるための代表N値について

代表 N 値（設計用採用値）の決め方 ― 要点まとめ 1．実務でよく尋ねられる「出典」の悩み（改訂） 発注者から「この整理手法の出典はどこか」と問われることが多い一方で、道路土工・道路橋示方書・地盤調査法など主要文献には 「平均値±1/2σで採用する」と明記した箇所はほとんど見当たりません。 しかし、 『報告書作成マニュアル［土質編］』（社団法人 全国地質調査業協会連合会） には、 ばらつきデータから設計用値を定める 三手法 （ a： 平均、 b： 最小／最大、 c： 平均と標準偏差による補正）が明記されています。 特に c) として次式の例が示され、 設計値 ＝ 平均値 − (1/2) × 標準偏差 （平均値を標準偏差で補正して設計値とする） と記載されています（第2章 2.7.3、p.74）。 したがって、本稿で用いる 「平均 − 1/2σ」 という整理は単なる経験則ではなく、 公的マニュアルに裏づけのある実務手法 です。 なお同マニュアルは、標準偏差の重み係数（1/2 など）の選び方には議論があること、すなわち k の選定はデータ量・ばらつき・重要度等に応じて調整する余地があることも併記しています。 出典： 『報告書作成マニュアル［土質編］』第2章 2.7.3「設計用土質定数の設定方法」p.74（a〜c の三手法と「平均 − (1/2) 標準偏差」の例示） 2．建築基礎構造設計規準・同解説（1974）の記述 地盤の場合、土質試験などによって得られる諸数値にかなりのばらつきがあるのがむしろ普通であり、このようにばらつきのある測定値の取り扱いについては、ばらつきの原因とその程度、基礎底面と各試料採取位置との相対的な関係など、複雑な要素が絡むため一概には言えない。しかし通常、1地点で同一土質層について得られた測定値は、その平均値を用いて差し支えないとする。そのうえで、とくにばらつきが著しい場合などは、平均値よりさらに安全側を見込んで 設計用採用値 ＝ （平均値） ± （標準偏差）／2 とする。ここで複号 ± は、算定結果が安全側となる符号をとる...
続きを読む