[Blender3.3] 日本の地形モデルメッシュ追加アドオンの開発 (6): 取り敢えずの完成

前回、スクリプトのファイルを分割して、それぞれの機能が見通し良くなりました。

そこで最後の仕上げに入ります。

まず、距離と緯度との変換のためにヒュベニの公式のモジュールを作ります。ファイル名は、hubeny_formyla.py です。

# -*- coding: utf-8 -*-
#
# package: terrain_model
# file: hubeny_formula.py

# Hubeny's Formula
#
# D = sqrt(Dy * M)^2 + (Dx * N * cos(P))^2)
# D: distance (m)
# Dy: diff of latitude (radian)
# Dx: diff of longitude (radian)
# P: average of latitude (radian)
# M = Rx * (1 - E^2) / W^3
# W = sqrt(1 - E^2 * sin(P)^2)
# N = Rx / W
# E = sqrt((Rx^2 - Ry^2) / Rx^2)
#
#       long radius(Rx), short radius(Ry)
# WGS84: 6,378,137.000, 6,356,752.314 245
#

import math

Rx = 6378137.0
Ry = 6356752.314245
E2 = (Rx ** 2 - Ry ** 2) / Rx ** 2 
#W = math.sqrt(1 - E2 * math.sin(P) ** 2)
#N = Rx / W

def angle_from(distance, latitude):
    '''
    parameter:
        distance (m)
        latitude (radian)
    '''
    W = math.sqrt(1 - E2 * math.sin(latitude) ** 2)
    N = Rx / W
    return distance / N / math.cos(latitude)

def distance_from(angle, latitude):
    '''
    parameter:
        angle (radian)
        latitude (radian)
    '''
    W = math.sqrt(1 - E2 * math.sin(latitude) ** 2)
    N = Rx / W
    return angle * N * math.cos(latitude)

そして、これを使ってモデル作成の処理を見直します。今までは、指定された緯度経度に該当するタイルを特定し、そのタイルの情報だけで地形モデルを作成していましたが、当初構想通り、指定された緯度経度に対し、指定された範囲にある地図タイルにピクセル単位でアクセスして、標高データを集めてメッシュを作成します。

そのように改修した model.py を示します。

# -*- coding: utf-8 -*-
#
# package: terrain_model
# file: model.py

import bpy
import math
import requests
from bpy_extras.object_utils import object_data_add
from mathutils import Vector
from . import hubeny_formula
from trace import Trace

def T(str):
    return bpy.app.translations.pgettext(str)

def add_object(self, context):
    
    lat = self.prop_lat
    lon = self.prop_lon
    zoom = self.prop_zoom
    scope = self.prop_scope # (km)
    Trace.print("add_object(lat=%f, lon=%f, scope=%f, zoom=%d)" % (lat, lon, scope, zoom))
    
    scope_x = hubeny_formula.angle_from(scope * 1000, lat) / (2 * math.pi)

    (x, y) = TerrainTile.latlon_to_xy(lat, lon)
    center_px = x * (256 << zoom)
    center_py = y * (256 << zoom)
    scope_p = scope_x * (256 << zoom)
    start_px = math.ceil(center_px - scope_p)
    start_py = math.ceil(center_py - scope_p)
    end_px = math.floor(center_px + scope_p)
    end_py = math.floor(center_py + scope_p)
    tile = get_terrain_tile(start_px, start_py, end_px, end_py, zoom)

    scene = context.scene
    scale_length = scene.unit_settings.scale_length

    start_x = hubeny_formula.distance_from((start_px - center_px) / (256 << zoom) * 2 * math.pi, lat) / scale_length
    start_y = hubeny_formula.distance_from((center_py - start_py) / (256 << zoom) * 2 * math.pi, lat) / scale_length
    width = hubeny_formula.distance_from((end_px - start_px) / (256 << zoom) * 2 * math.pi, lat) / scale_length
    height = hubeny_formula.distance_from((start_py - end_py) / (256 << zoom) * 2 * math.pi, lat) / scale_length
    mesh = tile.create_mesh(start_x, start_y, width, height)

    # useful for development when the mesh may be invalid.
    # mesh.validate(verbose=True)
    object_data_add(context, mesh, operator = self)

tile_map = {}

def get_terrain_tile(start_px: int, start_py: int, end_px: int, end_py: int, zoom: int):
    alts = []
    for py in range(start_py, end_py):
        ty = py >> 8
        tpy = py & 255
        for px in range(start_px, end_px):
            tx = px >> 8
            tpx = px & 255
            s = (zoom, tx, ty)
            if s in tile_map.keys():
                tile = tile_map[s]
            else:
                tile = get_cyberjapandata(zoom, tx, ty)
                tile_map[s] = tile
            alts.append(tile.altitude_at(tpx, tpy))
    
    tile = TerrainTile(pixel_x = end_px - start_px, pixel_y = end_py - start_py, altitudes = alts)
    return tile

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, pixel_x = PIXEL_NUM, pixel_y = PIXEL_NUM, altitudes = [], csvtext = ""):
        self.altitudes = altitudes
        self.pixel_x = pixel_x
        self.pixel_y = pixel_y
        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_y = self.pixel_y
        pixel_x = self.pixel_x
        lines = csvtext.split('\n')
        if len(lines) < pixel_y:
            raise ValueError("num of lines is less than pixel_y(%d): %d lines" % (pixel_y, len(lines)))
        
        for i, line in enumerate(lines[:pixel_y]):
            cols = line.split(',')
            if len(cols) != pixel_x:
                raise ValueError("num of columns at line %d is not pixel_x(%d): %d cols" % (i + 1, pixel_x, len(cols)))
            for col in cols:
                col = col.strip()
                if col == "e":
                    alts.append(None)
                else:
                    alts.append(float(col))
        self.altitudes = alts
        
    def altitude_at(self, x, y):
        if (x < 0 or x >= self.pixel_x):
            return None
        if (y < 0 or y >= self.pixel_y):
            return None
        index = y * self.pixel_x + x
        if (len(self.altitudes) < index):
            return None
        return self.altitudes[index]
        
    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_y):
            y = start_y + (index_y + 0.5) * height / self.pixel_y
            
            for index_x, alt in enumerate(self.altitudes[(index_y * self.pixel_x):((index_y + 1) * self.pixel_x)]):
                if alt is not None:
                    x = start_x + (index_x + 0.5) * width / self.pixel_x
                    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_y - 1):
            iy0 = index_y * self.pixel_x
            iy1 = iy0 + self.pixel_x
            for ix0, ix1 in enumerate(range(1, self.pixel_x)):
                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

機能面での改修は以上ですが、新しくhubeny_formula.pyを作成したので、これのreloadをするように、__init__.py の __all__ に、hubeny_formula を追加します。

# -*- coding: utf-8 -*-
#
# package: terrain_model
# file: __init__.py

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",
}

__all__ = [ 'translation', 'operator', 'model', 'hubeny_formula' ]

from . import translation
from . import operator

import bpy
from trace import Trace

## 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.app.translations.register(__name__, translation.translation_dict)
    operator.register()
#    bpy.utils.register_manual_map(add_object_manual_map)


def unregister():
    Trace.print("unregister()")
    operator.unregister()
    bpy.app.translations.unregister(__name__)
#    bpy.utils.unregister_manual_map(add_object_manual_map)


if __name__ == "__main__":
    register()

これで、当初考えていた内容が実装できました。

実際に緯度経度を指定して拾ったデータでメッシュオブジェクトを作成した例を示します。

ぱっと見でお分かりになるでしょうか。富士山です。色が無いと味気ないですね。

色の付け方については、これから勉強していきたいと思います。

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