[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)】

このブログの人気の投稿

[Blender3.3] mmd_toolsはどれが最新？

 キャラクターのモデルを探すとMMD形式のものが色々とあります。 このMMD形式のファイルは、素のBlenderでは読み込めません。読む込むにはmmd_toolsが必要です。 で、調べているとmmd_toolsには複数あるようで、どれが良いのか混乱するので、調べた結果をまとめておきます。なお、今後も変化する可能性があるので、あくまで本日(2022/9/24)時点の情報です。 ・https://github.com/sugiany/blender_mmd_tools これが大元のようですが、7年くらい前で更新が停止。 ・https://github.com/powroupi/blender_mmd_tools 上記より新しいと思われるのがこれですが、これもまた更新が停止した模様です。これのホームでは、下記に移行したとのアナウンスが記載されていました。 ・https://github.com/UuuNyaa/blender_mmd_tools これが最新で有効なmmd_toolsみたいです。Blender3.3LTSに対応しています。インストール方法など、wikiにて詳細に述べられています。 MMD Tools | MMD & Blender Wiki | Fandom (このWikiページの最上部の表示は広告であり、記事の一部ではないので、十分ご注意を) また、読み込んだMMDデータの扱いについては、それぞれのデータの著作権表示その他を十分読んで理解した上で、それに従ってください。大体は個人的に研究する分には問題ないと思いますが、再配布や商用利用については制限がついていると思いますので、注意が必要です。
続きを読む

パズドラ 転生龍刀士降臨周回編成 (ベア戦)

イメージ
龍刀士のイナやセリカを転生させたいと思ったことはありませんか？ 思ったところで、これまで転生龍刀士降臨は、ハードルが高くて近寄れていませんでした。   が、正月ノルディスが手に入ったし、できるかなと編成を組んでみました。それが以下です。 (文章の続きは画像の後)   アシストも潜在覚醒も無効で、ペア戦なので超覚醒もなく、画面を見たまんまです。 下側の正月ノルディスの隣は、ライゼクス装備。なぜ装備を単独で置いているかと言えば、ダメージ吸収無効に属性吸収無効スキルを持ち、コンボドロップがある光属性の物がこれしかなかったからです。   コンボドロップはB2の夜想詠神ツクヨミ対策ですが、これでも不足。その解説の前に立ち回りを。 正月ノルディスは単にノルディスって書いてます。Lがリーダー、Hが助っ人、Sがサブにいるノルディスです。 B1: [上] Sノルディス、御坂美琴、Hノルディス、アーニャ達 (ずらし) B2: [下] シェザル (追加5コンボ) B3: [上] Sノルディス (ずらし) B4: [下] Lノルディス (ずらし) B5: [上] ゼウスGIGA (ずらし)   → [下] アーニャ達 (ずらし) B6: [上] 御坂美琴、サノス (ずらし) B7: [下] キオまたはアキネなら ライゼクス装備、ノルディス×3 (ずらし)            セリカなら 1コンボ×3、Hノルディス (ずらし)            キョウリまたはイナなら ノルディス×3 (ずらし) 面倒なのはB2で。コンボドロップが3つ残って巻き込めるなら追加2コンボですが、コンボドロップが落ちコンで減ってしまっていると、その分、コンボを組まねばなりません。まあ、7x6盤面なので、最悪丸々5コンボ組まなければならなくても、何とかなりそうではあります。   これで、一応龍刀士を転生進化させられるようにはなりました。 実際に転生させるのは必要性が生じてからかな、と思ってます。   ではでは。
続きを読む

パズドラ HUNTER☓HUNTER キメラアント周回編成 (ペア戦)

イメージ
いや、なんでしょう、これは？ 事件？いや、事故？ 最初は軽い気持ちで、あーコラボ始まったんだぁ〜って思ったんです。 それでいつものように、クリア報酬でコラボガチャ1回が付いてくるダンジョンをやろうかな、とまずはHUNTER×HUNTER 幻影旅団の上級をクリアし、次はキメラアントでもやるかと攻略サイトで調べ始めたんです。   すると、何故かサレーネをリーダーにしたランク上げ編成と題した物が紹介されてたんです。  コラボダンジョンでなんでサレーネ使うの？って思ったわけですよ。   うーん、分からんと思って、普通に編成組んでダンジョンに挑んでクリアして。成功報酬の画面でランクアップしたなーと思いました。 そう、まだ、この時点では異変に気付いていなかったんです。 で、記事書くために安定してクリアできるか確認するべかと何度かやったら、またランクアップした。   「ん？」   おかしい。そんなに簡単にランクアップするわけないのに。 で、良く見ると取得経験値の値が見慣れない値であることに気付きました。   25,000,000   にせんごひゃくまん？   たまに配信される経験値ボーナスダンジョンも取得経験値は50万ではなかったろうか。   うむ、なるほど、サレーネをリーダーにしようと考える訳だ。   ようやく異常事態に気付いて、サレーネをリーダーにした周回編成を検討したのでした。 それが下の物です。   (文章の続きは絵の下です) 潜在覚醒ですが、正月ノルディスにはダメージ上限解放を付けて、4連ガチャドラにはスキルブースト++を付けて、正月ノルディスのスキルレベルを1とか2で済むようにしました。   (下記ノルディスは、すべて正月ノルディス) B1:[上] 助っ人ノルディス、右ノルディス、ロー、サレーネ (ずらし) B2:[下] 右ノルディス(助っ人の左側)、助っ人ノルディス (ずらし※) B3:[上] 中央ノルディス、左ノルディス、助っ人ノルディス (ずらし) B4:[下] 4連ガチャドラ、虹村億泰&ザ・ハンド、左ノルディス、助っ人ノルディス (ずらし) ※: 安定さが微妙なのはB2の8コンボのところです。木と光で2コンボ、コンボドロップ生成スキルが4つで6コンボなので、あと3コンボ必要です。B1で落ちたコンボドロップ
続きを読む