[Blender3.3] 日本の地形モデルメッシュ追加アドオンの開発 (4): 地理院タイル(標高タイル)の読み込み

 さて、いよいよ地理院タイルを読み込んで、地形モデルを作成したいと思います。

地理院タイルタイルの仕様はこちら(標高タイルの作成方法と地理院地図で表示される標高値について)を参照してください。

とは言っても、いきなり仕様通りのメッシュを作るのはハードルが高いので、まずは入力された緯度経度とズーム値で特定される一つの標高タイルを取得してそれをメッシュ化することにします。

スクリプトのコードは以下のようにしました。

# -*- coding: utf-8 -*-
 
bl_info = {
    "name": "New Terrain Model (JP)",
    "author": "Shiki Kuraga",
    "version": (1, 0),
    "blender": (3, 3, 0),
    "location": "View3D > Add > Mesh > New Terrain Model (JP)",
    "description": "Adds a new Terrain Model (JP) Mesh Object",
    "warning": "",
    "support": "TESTING",
    "doc_url": "",
    "category": "Add Mesh",
}
 
translation_dict = {
    "ja_JP": {
        ("*", "Latitude"):
            "緯度",
        ("*", "Longitude"):
            "経度",
        ("*", "Scope (km)"):
            "範囲長さ (km)",
        ("*", "Terrain Model (JP)"):
            "地形モデル (日本)",
        ("*", "Add Terrain Model Mesh (JP)"):
            "地形モデルの追加 (日本)"
    }
}
 
import bpy
import math
import requests
from bpy.types import Operator
from bpy.props import FloatProperty, IntProperty
from bpy_extras.object_utils import AddObjectHelper, object_data_add
from mathutils import Vector
from trace import Trace
 
def T(str):
    return bpy.app.translations.pgettext(str)
 
def add_object(self, context):
    Trace.print("add_object()")
    
    scene = context.scene
    scale_length = scene.unit_settings.scale_length
 
    scale_x = self.prop_scope * 1000 / scale_length #self.scale.x
    scale_y = self.prop_scope * 1000 / scale_length #self.scale.y
 
    zoom = self.prop_zoom
    (x, y) = TerrainTile.latlon_to_xy(self.prop_lat, self.prop_lon)
    ix = int(x * (1 << zoom))
    iy = int(y * (1 << zoom))
    
    tile = get_cyberjapandata(zoom, ix, iy)
    mesh = tile.create_mesh(-1 * scale_x, 1 * scale_y, 2 * scale_x, -2 * scale_y)
    
#    verts = [
#        Vector((-1 * scale_x, 1 * scale_y, 0)),
#        Vector((1 * scale_x, 1 * scale_y, 0)),
#        Vector((1 * scale_x, -1 * scale_y, 0)),
#        Vector((-1 * scale_x, -1 * scale_y, 0)),
#    ]
#
#    edges = []
#    faces = [[0, 1, 2, 3]]
#
#    mesh = bpy.data.meshes.new(name = "Terrain Model Mesh")
#    mesh.from_pydata(verts, edges, faces)
    # useful for development when the mesh may be invalid.
    # mesh.validate(verbose=True)
    object_data_add(context, mesh, operator = self)
 
 
def get_cyberjapandata(zoom: int, x: int, y: int, dataid = "DEM10B"):
    """Get cyberjapandata"""
# URL：https://cyberjapandata.gsi.go.jp/xyz/dem5a/{z}/{x}/{y}.txt（DEM5A テキスト形式）zoom: 1...15
# URL：https://cyberjapandata.gsi.go.jp/xyz/dem5b/{z}/{x}/{y}.txt（DEM5B テキスト形式）zoom: 1...15
# URL：https://cyberjapandata.gsi.go.jp/xyz/dem/{z}/{x}/{y}.txt（DEM10B テキスト形式） zoom: 1...14
 
    # check dataid
    if dataid == "DEM10B":
        (dtype, zmax) = ("dem", 14)
    elif dataid == "DEM5A":
        (dtype, zmax) = ("dem5a", 15)
    elif dataid == "DEM5B":
        (dtype, zmax) = ("dem5b", 15)
    else:
        raise ValueError("Invalid DATAID: %s" % dataid)
        
    # check zoom
    if type(zoom) is not int:
        raise TypeError("param zoom isn't int: %s" % type(zoom))
    elif zoom < 1 or zoom > zmax:
        raise ValueError("Invalid value of z(1-%d): %d" % (zmax, zoom))
        
    # check x
    tilenum = (1 << zoom)
    if type(x) is not int:
        raise TypeError("param x isn't int: %s" % type(x))
    elif x < 0 or x >= tilenum:
        Trace.print("Out of range value x(0-%d): %d" % (tilenum, x))
        return TerrainTile()
 
    # check y
    if type(y) is not int:
        raise TypeError("param y isn't int: %s" % type(y))
    elif y < 0 or y >= tilenum:
        Trace.print("Out of range value y(0-%d): %d" % (tilenum, y))
        return TerrainTile()
    
    # get data
    url = "https://cyberjapandata.gsi.go.jp/xyz/%s/%d/%d/%d.txt" % (dtype, zoom, x, y)
    Trace.print("get data from url: <%s>" % url)
    res = requests.get(url)
    s = res.content.decode()
    tile = TerrainTile(csvtext = s)
    Trace.print(tile)
    return tile
 
 
class TerrainTile:
    """Data Handler of Terrain Tile"""
    
    PIXEL_NUM = 256
    
    def latlon_to_xy(latitude,longitude):
        x = ((math.pi + longitude) / (2 * math.pi)) % 1
        
        siny = math.sin(latitude)
        siny = max(-0.999, min(0.999, siny))
        y = 0.5 - math.log((1 + siny) / (1 - siny)) / (4 * math.pi)
        
        return (x, y)
    
    def __init__(self, csvtext = ""):
        self.altitudes = []
        if len(csvtext) > 0:
            self.load_from_csvtext(csvtext)
 
    def __str__(self):
        return "" % len(self.altitudes)
        
    def load_from_csvtext(self, csvtext):
        alts = []
        if len(csvtext.strip()) == 0:
            self.altitudes = alts
            return
        PIXEL_NUM = TerrainTile.PIXEL_NUM
        lines = csvtext.split('\n')
        if len(lines) < PIXEL_NUM:
            raise ValueError("num of lines is less than PIXEL_NUM(%d): %d lines" % (PIXEL_NUM, len(lines)))
        
        for i, line in enumerate(lines[:PIXEL_NUM]):
            cols = line.split(',')
            if len(cols) != PIXEL_NUM:
                raise ValueError("num of columns at line %d is not PIXEL_NUM(%d): %d cols" % (i + 1, PIXEL_NUM, len(cols)))
            for col in cols:
                col = col.strip()
                if col == "e":
                    alts.append(None)
                else:
                    alts.append(float(col))
        self.altitudes = alts
        
    def create_mesh(self, start_x, start_y, width, height, name = "Terrain Model Mesh", adjust = True):
        offset_z = 0
        if adjust:
            max_z = max(filter(None, self.altitudes))
            min_z = min(filter(None, self.altitudes))
            Trace.print("min_z = ", min_z)
            if min_z > 0:
                offset_z = -min_z
            elif max_z < 0:
                offset_z = -max_z
            
        idmap = []
        verts = []
        for index_y in range(0, self.PIXEL_NUM):
            y = start_y + (index_y + 0.5) * height / self.PIXEL_NUM
            
            for index_x, alt in enumerate(self.altitudes[(index_y * self.PIXEL_NUM):((index_y + 1) * self.PIXEL_NUM)]):
                if alt is not None:
                    x = start_x + (index_x + 0.5) * width / self.PIXEL_NUM
                    idmap.append(len(verts))
                    verts.append(Vector((x, y, alt + offset_z)))
                else:
                    idmap.append(-1)
        faces = []
        for index_y in range(0, self.PIXEL_NUM - 1):
            iy0 = index_y * self.PIXEL_NUM
            iy1 = iy0 + self.PIXEL_NUM
            for ix0, ix1 in enumerate(range(1, self.PIXEL_NUM)):
                face = [idmap[i] for i in [ix0 + iy0, ix1 + iy0, ix1 + iy1, ix0 + iy1] if idmap[i] >= 0]
                if (len(face) >= 3):
                    faces.append(face)
        
        if len(faces) > 0:
            edges = []
            mesh = bpy.data.meshes.new(name = name)
            mesh.from_pydata(verts, edges, faces)
            return mesh
        
        return None
class OBJECT_OT_add_object2(Operator, AddObjectHelper):
    """Create a new Mesh Object2"""
    bl_idname = "mesh.add_object2"
    bl_label = T("Add Terrain Model Mesh (JP)")
    bl_options = {'REGISTER', 'UNDO'}
    prop_lat: FloatProperty(
        name = "Latitude",
        description = "Latitude",
        default = 35.0 * math.pi / 180.0,
        subtype = 'ANGLE',
        min = -0.495 * math.pi,
        max = 0.495 * math.pi
    )
    prop_lon: FloatProperty(
        name = "Longitude",
        description = "Longitude",
        default = 135.0 * math.pi / 180.0,
        subtype = 'ANGLE',
        min = -1 * math.pi,
        max = 1 * math.pi
    )
    prop_scope: FloatProperty(
        name = "Scope (km)",
        description = "Scope (km)",
        default = 1.0,
        step = 10, # actual value = 10/100
        min = 0.1
    )
    prop_zoom: IntProperty(
        name = "Zoom",
        description = "Zoom",
        default = 14,
        min = 3,
        max = 15
    )
    
    def __init__(self):
        Trace.print("__init__")
        
    def __del__(self):
        Trace.print("__del__")
    def execute(self, context):
        add_object(self, context)
        return {'FINISHED'}
    def invoke(self, context, event):
        
        wm = context.window_manager
        
        # invoke properties dialog for this object
        return wm.invoke_props_dialog(self)
# Registration
def add_object_button2(self, context):
    # insert separator
    self.layout.separator()
    self.layout.operator(
        OBJECT_OT_add_object2.bl_idname,
        text = T("Terrain Model (JP)"),
        icon = 'PLUGIN')
## This allows you to right click on a button and link to documentation
#def add_object_manual_map():
#    url_manual_prefix = "https://docs.blender.org/manual/en/latest/"
#    url_manual_mapping = (
#        ("bpy.ops.mesh.add_object2", "scene_layout/object/types.html"),
#    )
#    return url_manual_prefix, url_manual_mapping
def register():
    Trace.print("register()")
    bpy.utils.register_class(OBJECT_OT_add_object2)
    bpy.app.translations.register(__name__, translation_dict)
#    bpy.utils.register_manual_map(add_object_manual_map)
    bpy.types.VIEW3D_MT_mesh_add.append(add_object_button2)
def unregister():
    Trace.print("unregister()")
    bpy.app.translations.unregister(__name__)
    bpy.utils.unregister_class(OBJECT_OT_add_object2)
#    bpy.utils.unregister_manual_map(add_object_manual_map)
    bpy.types.VIEW3D_MT_mesh_add.remove(add_object_button2)
if __name__ == "__main__":
    # to avoid duplicate registration error of translation_dict
    if (Trace.activated):
        bpy.app.translations.unregister(__name__)
        
    Trace.on()
    register()

ソースコードを超簡単に解説しますと、TerrainTileというタイル情報を保持するクラスを定義しています。

標高タイルを取得する関数が get_cyberjapandata() で、これの結果がTerrainTileで戻るとしてあります。

TerrainTileには、自分で保持しているデータをもとにメッシュを作成するcreate_mesh()という関数を持っているので、これでメッシュに変換して追加するようにしたということです。

TerrainTileを作ったのは、ポップアップ画面で指定した領域のデータを入れるのにも使えるようにすれば、create_mesh()の処理を共通的に使えて楽かなと考えたためです。

こうして作成したスクリプトを実行し、北緯35度、東経135度を指定して標高タイルを取得した結果を下に示します。

それらしく読み込めたように見えます。

これはDEM10Bなので、ピクセル一つで約10m、タイル一つは約2.5km四方でしょうか。本来ポップアップ画面では範囲1kmを指定したので2km四方だから気持ち狭い範囲となるわけですが、大体こんな感じになるのだろいうイメージが持てました。

次は仕様通りの読み込み範囲にすることをと思いつつ、開発環境の整備をしないとかな、と考え中です。

【前回(3)】日本の地形モデルメッシュ追加アドオンの開発 (4) 【次回(5)】