【QGIS】国土地理院ベクターのズーム制限解除 ＆ 巨大点群からの傾斜区分図自動作成（岐阜県バージョン）

【QGIS】国土地理院ベクターのズーム制限解除 ＆ 巨大点群からの傾斜区分図自動作成（岐阜県バージョン）

1. 国土地理院ベクタータイルの「拡大すると消える」問題を解決する

QGISで「国土地理院ベクター（最適化ベクトルタイル）」を表示した際、標準スタイルではズームレベル16.1以上に拡大すると建物や道路が消えてしまう仕様になっています。これでは詳細な図面作成ができません。

ベクタタイルレイヤはプロパティから一括でズームレベルを変更できませんが、Pythonコンソールを使えば一瞬で全項目の制限を解除できます。

【解決手順】

  1. レイヤパネルで「国土地理院ベクター」レイヤを選択。

  2. メニューの [プラグイン] ＞ [Pythonコンソール] を開く。

  3. 下部の入力欄に以下の1行を貼り付けてEnter。

Python
layer = iface.activeLayer(); styles = layer.renderer().styles(); [s.setMaxZoomLevel(22) for s in styles]; layer.renderer().setStyles(styles); layer.triggerRepaint()

これでレベル22まで拡大しても地物が表示され続けます。

2. 巨大点群から傾斜区分図を一括作成する（岐阜県仕様）

数千万点のデータから傾斜区分図を作る際、QGISのGUI操作ではフリーズのリスクが高まります。ここでは、用途に応じて使い分けられる「超高速版」と「高精度フルバージョン」の2つのスクリプトを紹介します。

傾斜区分の定義（5段階）

  • 値 1：0° - 2°（緑色）

  • 値 2：2° - 15°（黄緑色）

  • 値 3：15° - 30°（薄い黄色）

  • 値 4：30° - 45°（山吹色）

  • 値 5：45° - 90°（オレンジ色）

A. 【超高速版】とりあえず全体を素早く確認したい時

1000万点超えのデータでも数分で完了します。微地形は多少簡略化されますが、スピード重視の現場確認に最適です。

Python
import os
import numpy as np
import tkinter as tk
from tkinter import filedialog
try:
    from osgeo import gdal, ogr, osr
except ImportError:
    print("エラー: GDALライブラリが見つかりません。OSGeo4W Shellから実行してください。")
    import sys
    sys.exit()

def process_fast_slope():
    root = tk.Tk()
    root.withdraw()
    root.attributes('-topmost', True)
    input_path = filedialog.askopenfilename(title="処理するファイルを選択", filetypes=[("GeoPackage/Shapefile", "*.gpkg;*.shp")])
    if not input_path: return

    base_dir = os.path.dirname(input_path)
    base_name = os.path.splitext(os.path.basename(input_path))[0]
    output_name = base_name.replace("_clipped", "")

    dem_tif = os.path.join(base_dir, f"{output_name}_dem.tif")
    slope_tif = os.path.join(base_dir, f"{output_name}_slope.tif")
    reclass_tif = os.path.join(base_dir, f"{output_name}_reclass.tif")
    
    z_field = "Z"      
    res = 1.0          

    ds = ogr.Open(input_path)
    layer = ds.GetLayer()
    x_min, x_max, y_min, y_max = layer.GetExtent()
    srs = layer.GetSpatialRef()
    srs_wkt = srs.ExportToWkt() if srs else 'EPSG:6675' 
    
    print(f"超高速処理開始: {os.path.basename(input_path)}")
    
    # 1. ラスタライズ
    gdal.Rasterize(dem_tif, input_path, options=gdal.RasterizeOptions(
        format='GTiff', outputBounds=[x_min, y_min, x_max, y_max],
        xRes=res, yRes=res, attribute=z_field, noData=-9999, outputType=gdal.GDT_Float32, outputSRS=srs_wkt))

    # 2. 隙間補間
    ds_dem = gdal.Open(dem_tif, gdal.GA_Update)
    gdal.FillNodata(targetBand=ds_dem.GetRasterBand(1), maskBand=None, maxSearchDist=20, smoothingIterations=0)
    ds_dem = None

    # 3. 傾斜計算
    gdal.DEMProcessing(slope_tif, dem_tif, 'slope', options=gdal.DEMProcessingOptions(format='GTiff', computeEdges=True, slopeFormat='degree'))

    # 4. 再分類
    ds_slope = gdal.Open(slope_tif)
    band = ds_slope.GetRasterBand(1)
    arr = band.ReadAsArray()
    nodata = band.GetNoDataValue()
    
    reclass = np.full(arr.shape, -9999, dtype=np.int16)
    v = (arr != nodata)
    reclass[v & (arr >= 0) & (arr <= 2)] = 1
    reclass[v & (arr > 2) & (arr <= 15)] = 2
    reclass[v & (arr > 15) & (arr <= 30)] = 3
    reclass[v & (arr > 30) & (arr <= 45)] = 4
    reclass[v & (arr > 45)] = 5

    driver = gdal.GetDriverByName('GTiff')
    ds_out = driver.Create(reclass_tif, ds_slope.RasterXSize, ds_slope.RasterYSize, 1, gdal.GDT_Int16)
    ds_out.SetGeoTransform(ds_slope.GetGeoTransform())
    ds_out.SetProjection(ds_slope.GetProjection())
    ds_out.GetRasterBand(1).WriteArray(reclass)
    ds_out.GetRasterBand(1).SetNoDataValue(-9999)
    ds_out = ds_slope = None
    print(f"✅ 完了: {reclass_tif}")

if __name__ == "__main__":
    process_fast_slope()
input("Enterキーを押して終了...")

B. 【高精度フルバージョン】公式図面・報告書作成用

全ての点を精密に結ぶ「線形補間」を行います。1000万点超えでは10時間以上かかる場合がありますが、微地形まで忠実に再現します。退勤前の実行を推奨します。

Python
import os
import numpy as np
import tkinter as tk
from tkinter import filedialog
try:
    from osgeo import gdal, ogr, osr
except ImportError:
    print("エラー: GDALライブラリが見つかりません。")
    import sys
    sys.exit()

def process_full_slope():
    root = tk.Tk()
    root.withdraw()
    root.attributes('-topmost', True)
    input_path = filedialog.askopenfilename(title="処理するファイルを選択", filetypes=[("GeoPackage/Shapefile", "*.gpkg;*.shp")])
    if not input_path: return

    base_dir = os.path.dirname(input_path)
    base_name = os.path.splitext(os.path.basename(input_path))[0]
    output_name = base_name.replace("_clipped", "")

    dem_tif = os.path.join(base_dir, f"{output_name}_dem.tif")
    slope_tif = os.path.join(base_dir, f"{output_name}_slope.tif")
    reclass_tif = os.path.join(base_dir, f"{output_name}_reclass.tif")
    
    z_field = "Z"      
    res = 1.0          

    ds = ogr.Open(input_path)
    layer = ds.GetLayer()
    x_min, x_max, y_min, y_max = layer.GetExtent()
    width = int((x_max - x_min) / res)
    height = int((y_max - y_min) / res)
    srs = layer.GetSpatialRef()
    srs_wkt = srs.ExportToWkt() if srs else 'EPSG:6675'
    
    print(f"高精度版処理開始: {os.path.basename(input_path)}")
    
    # 1. 線形補間(Grid)
    gdal.Grid(dem_tif, input_path, options=gdal.GridOptions(
        format='GTiff', outputBounds=[x_min, y_min, x_max, y_max],
        width=width, height=height, algorithm='linear', zfield=z_field, outputSRS=srs_wkt))

    # 2. 傾斜計算
    gdal.DEMProcessing(slope_tif, dem_tif, 'slope', options=gdal.DEMProcessingOptions(format='GTiff', computeEdges=True, slopeFormat='degree'))

    # 3. 再分類
    ds_slope = gdal.Open(slope_tif)
    band = ds_slope.GetRasterBand(1)
    arr = band.ReadAsArray()
    nodata = band.GetNoDataValue()
    
    reclass = np.full(arr.shape, -9999, dtype=np.int16)
    v = (arr != nodata)
    reclass[v & (arr >= 0) & (arr <= 2)] = 1
    reclass[v & (arr > 2) & (arr <= 15)] = 2
    reclass[v & (arr > 15) & (arr <= 30)] = 3
    reclass[v & (arr > 30) & (arr <= 45)] = 4
    reclass[v & (arr > 45)] = 5

    driver = gdal.GetDriverByName('GTiff')
    ds_out = driver.Create(reclass_tif, ds_slope.RasterXSize, ds_slope.RasterYSize, 1, gdal.GDT_Int16)
    ds_out.SetGeoTransform(ds_slope.GetGeoTransform())
    ds_out.SetProjection(ds_slope.GetProjection())
    ds_out.GetRasterBand(1).WriteArray(reclass)
    ds_out.GetRasterBand(1).SetNoDataValue(-9999)
    ds_out = ds_slope = None
    print(f"✅ 完了: {reclass_tif}")

if __name__ == "__main__":
    process_full_slope()
input("Enterキーを押して終了...")

おわりに

作成された _reclass.tif をQGISに読み込み、「単バンド疑似カラー」「正確」モードで5色を割り当てれば完成です。

コメント

  1. eon2026年5月7日 18:55

    建物の色を砂防基盤図に合わせる。
    QColor(170, 5, 134)

    from qgis.PyQt.QtGui import QColor
    from qgis.core import Qgis

    layer = iface.activeLayer()
    styles = layer.renderer().styles()
    pink = QColor(170, 5, 134)

    for s in styles:
    # ここを layerName() に修正しました
    if s.layerName() == 'building':
    sym = s.symbol()
    # 線（LineString）の場合は線の色を変更
    if sym.type() == Qgis.SymbolType.Line:
    sym.setColor(pink)
    # 面（Polygon）の場合は外枠（Stroke）の色を変更
    elif sym.type() == Qgis.SymbolType.Fill:
    sym.symbolLayer(0).setStrokeColor(pink)

    layer.renderer().setStyles(styles)
    layer.triggerRepaint()

    返信削除
  2. eon2026年5月7日 18:55

    ズームレベルを破壊（拡大すると、非表示になる）

    layer = iface.activeLayer(); styles = layer.renderer().styles(); [s.setMaxZoomLevel(22) for s in styles]; layer.renderer().setStyles(styles); layer.triggerRepaint()

    返信削除

コメントを投稿

このブログの人気の投稿

IJCAD（AutoCAD互換CAD） コマンド早見表

AutoCAD LT 2025 コマンド早見表（完全版） AutoCAD LT 2025 主要コマンド早見表【日本語備考付き・完全版】 AutoCAD LT で頻出するコマンドを アルファベット順 にまとめました。 備考欄には「実務で役立つワンポイント」を一行で追加。★印はブロック・外部参照まわりで特に押さえておきたいコマンドです。 略号 正式コマンド 機能概要 備考 A ARC 円弧を作成 半径・３点・開始終了角で指定可 AA AREA 面積と周長を計算 ポリライン選択で複数領域合算 ADC ADCENTER コンテンツセンター ブロックや画像をドラッグ挿入 AL ALIGN 2D/3D で位置合わせ 3 点指定でスケール同時調整 AR ARRAY 配列複写 楕円・パス配列も可能 AS ASSISTANTOPEN Autodesk Assistant Web ヘルプ検索に便利 ATI ATTIPEDIT★ ブロック属性文字を直接編集 図面上でダブルクリックと同等 ATT ATTDEF★ ブロック属性を定義 タグ・既定値・プロンプトを設定 ATE ATTEDIT★ 挿入済み属性を編集 EATTEDIT GUI と同等 B BLOCK★ ブロック定義を作成 書き出し前に 基点 を忘れず指定 BC BCLOSE★ ブロックエディタを閉じる SAVEAS 不可状態から復帰 BE BEDIT★ ブロックエディタを開く 動的ブロック対応の本格編集 BH HATCH ハッチング 「固有原点」でパターンズレ防止 BO BOUNDARY 領域/ポリライン生成 内部島検出で高速ポリライン化 BR BREAK オブジェクトを分割 F オプションで１点除去 BS BSAVE★ ブロック定義を保存 BEDIT 中の上書き保存 BVS BVSTATE★ 動的ブロック可視状態 仕上違いを１定義で管理 C CIRCLE 円を作成 2P/3P 補助で簡易作図 CBAR CONSTRAIN...
続きを読む

設計定数を求めるための代表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 とする。ここで複号 ± は、算定結果が安全側となる符号をとる...
続きを読む

日本でよく使う 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 福岡・佐賀・熊本・大分・宮崎、鹿児島...
続きを読む