日々適当

正距円筒図法の日本地図を他の図法に変形したいわけですけど、Pythonのモジュールとしてはcartopyとmatplotlibの組み合わせでできそうだと思って格闘したんですけど、荒れた絵しか出てこないのですね。単に無知な結果なのか、そもそも無理なのか…
なので諦めて、3DCGソフトのグリッドに日本地図を貼り付けて、それを変形させることにしました。

Pythonが使えればいいのでなんでもいいんですけど、せっかくなのでBlenderで。座標変換にはPyprojモジュールを使います。
macOS上で pip3 install pyproj でOSのPythonにはそのモジュールが入りましたが、BlenderのPythonでは認識されないので一手間加えます。

Blender上で
import sys
sys.path.append('/usr/local/lib/python3.7/site-packages')
をするだけなんですけどね。これでpyprojが使えるようになりました。

ということで、突貫工事の稚拙なスクリプト。1000x1000のグリッドで実行完了まで数秒かかりますw

import bpy
from mathutils import Vector, Matrix
from pyproj import Proj
import numpy as np

#変形後の形状の東西方向のサイズ(m)
gridWidth = 100.0

#日本各地の端っこ（Google mapでクリックして返ってくる値w）
north = [ 45.523015 , 141.936587 ] #宗谷岬
#east= [ 43.385437 , 145.817579 ] #納沙布岬
east = [ 45.509032 , 148.893014 ] #択捉島
south = [ 24.045030 , 123.791862 ] #波照間島
west = [ 24.449230 , 122.932030 ] #与那国島

#中心地としてとりあえず明石を指定してみる。
akashi = [ 34.0 + 39.0/60.0, 135.0  ] #明石

#ランベルト正角円錐図法
#標準緯線が何が正解か分からんので、とりあえず国土地理院のページの地図の一つを参考にしている。
#https://www.gsi.go.jp/chubu/minichishiki22.html
lat_1 = 34.0 + 5.0/60.0
lat_2 = 37.0 + 54.0/60.0 + 30.0/3600.0
converter = Proj( proj='lcc', R=6371200,  lat_1=lat_1, lat_2=lat_2, lon_0=akashi[1], ellps='WGS84' )

obj = bpy.context.object
mesh = obj.data
verts = mesh.vertices

xMax = -1000000.0
yMax = -1000000.0
xMin = 1000000.0
yMin = 1000000.0
vertsList = []

for vert in verts:
    vertPos =  vert.co[:2] 
    vertsList.append( list( vertPos ) )
    if vertPos[0] > xMax:
        xMax = vertPos[0]
    if vertPos[0]  < xMin:
        xMin = vertPos[0] 
    if vertPos[1] > yMax:
        yMax = vertPos[1]
    if vertPos[1] < yMin:
        yMin = vertPos[1]

newXMax = -10000000.0
newYMax = -10000000.0
newXMin = 10000000.0
newYMin = 10000000.0
i = 0

cntX , cntY  = converter( akashi[1], akashi[0] ) #中心地は明石にしている

for vertPos in vertsList:
    lon = ( vertPos[0] - xMin )/( xMax - xMin ) * ( east[1] - west[1] ) + west[1]
    lat = ( vertPos[1] - yMin )/( yMax - yMin ) * ( north[0] - south[0] ) + south[0]
    
    x , y = converter( lon, lat )
    x = x * 0.001
    y = y * 0.001
    
    y = y - cntY * 0.001
    
    if x > newXMax:
        newXMax = x
    if y > newYMax:
        newYMax = y
    if x < newXMin:
        newXMin = x
    if y < newYMin:
        newYMin = y
    
    v = Vector( ( x, y, 0.0 ) )
    verts[i].co = v
    i+=1

sclVal = gridWidth / (newXMax - newXMin)
obj.scale = ( sclVal, sclVal, sclVal )