2017-10-10

標高データに土地利用細分メッシュを重ねる（３）

ジオクラ

いろいろ調整して不自然さを除いてみた（農地を「草」ブロックにした。一部、そのまま）。

f:id:pycra:20171010014735p:plain

地面に立っての風景（本当はこの草じゃなくて普通の草を生やしたい）。 f:id:pycra:20171010014752p:plain

ちょっと高い位置から。奥の赤いブロックは「建物」がある場所。 f:id:pycra:20171010014805p:plain

スクリプト

# coding: utf-8

import mcpi.minecraft as minecraft
import mcpi.block as block
import numpy as np
import pycra.npydata as nd
import pycra.wood as wood
import random
import sys
from time import sleep

NPY_2ND_FILE = "npydata-{0:04d}-{1:02d}.npy"
NPY_3RD_FILE = "npydata-{0:04d}-{1:02d}-{2:02d}.npy"
L03B_2ND_FILE = "npyl03-b-{0:04d}-{1:02d}.npy"
L03B_3RD_FILE = "npyl03-b-{0:04d}-{1:02d}-{2:02d}.npy"

if len(sys.argv) < 3:
    sys.exit("Usage: loadnpy4.py mesh1 mesh2 [mesh3]")

mc = minecraft.Minecraft()

mesh1 = int(sys.argv[1])
mesh2 = int(sys.argv[2])
mesh3 = None
if len(sys.argv) >= 4:
    mesh3 = int(sys.argv[3])

if mesh3 is None:
    sys.exit("No implement !")

h_array = np.load(nd.NPY_DIR+"/"+NPY_3RD_FILE.format(mesh1, mesh2, mesh3))
l_array = np.load(nd.NPY_DIR+"/"+L03B_3RD_FILE.format(mesh1, mesh2, mesh3))

for z in range(150):
    for x in range(225):
        r = random.random()
        hval = h_array[z][x]
        lval = l_array[z][x]
        hval *= 0.2
        mc.setBlocks(x, hval-1, z, x, hval-7, z, block.DIRT)
        if lval == 100:  # 田
            mc.setBlock(x, hval, z, block.GRASS)
            if r < 0.1:
                mc.setBlock(x, hval+1, z, block.FERN)
            #mc.setBlock(x, hval, z, block.FARMLAND)
            #mc.setBlock(x, hval+1, z, block.WHEAT)
        elif lval == 200: # その他農地
            mc.setBlock(x, hval, z, block.GRASS)
            if r < 0.1:
                mc.setBlock(x, hval+1, z, block.FERN)
            #mc.setBlock(x, hval, z, block.HARDENED_CLAY_STAINED_ORANGE)
        elif lval == 500: # 森林
            mc.setBlock(x, hval, z, block.GRASS)
            if r < 0.03:
                wood.wood(mc, x, hval+1, z)
            elif r < 0.1:
                mc.setBlock(x, hval+1, z, block.FERN)
        elif lval == 600: # 荒地
            mc.setBlock(x, hval, z, block.HARDENED_CLAY_STAINED_YELLOW)
        elif lval == 700: # 建物
            mc.setBlock(x, hval, z, block.HARDENED_CLAY_STAINED_RED)
        elif lval == 901: # 道路
            mc.setBlock(x, hval, z, block.HARDENED_CLAY_STAINED_LIGHT_GRAY)
        elif lval == 902: # 鉄道
            mc.setBlock(x, hval, z, block.HARDENED_CLAY_STAINED_GRAY)
        elif lval == 1000: # その他用地
            mc.setBlock(x, hval, z, block.HARDENED_CLAY_STAINED_PINK)
        elif lval == 1100: # 河川及び湖沼
            mc.setBlock(x, hval, z, block.HARDENED_CLAY_STAINED_BLUE)
        elif lval == 1400: # 海浜
            mc.setBlock(x, hval, z, block.HARDENED_CLAY_STAINED_CYAN)
        elif lval == 1500: # 海水域
            mc.setBlock(x, hval, z, block.HARDENED_CLAY_STAINED_PURPLE)
        elif lval == 1600: # ゴルフ場
            mc.setBlock(x, hval, z, block.HARDENED_CLAY_STAINED_LIME)
        else:
            mc.setBlock(x, hval, z, block.DIRT)
    sleep(0.02)

2017-10-10

標高データに土地利用細分メッシュを重ねる（２）

ジオクラ

さすがに色分けするとあからさまなので、「森林」の場所に「木」を生やしてみる。

f:id:pycra:20171010012440p:plain

なんか、いい感じ。

「木」を減らして「草」を生やして、「田」の場所に「耕地」ブロックを配し、「小麦」ブロックを乗せてみた。

f:id:pycra:20171010012534p:plain

これはイマイチ。

木を減らしすぎたのと、「耕地」がハッキリしすぎて地形に合ってない。

やはり、この辺は標高データと合わせて調整したい。

スクリプト

# coding: utf-8

import mcpi.minecraft as minecraft
import mcpi.block as block
import numpy as np
import pycra.npydata as nd
import pycra.wood as wood
import random
import sys
from time import sleep

NPY_2ND_FILE = "npydata-{0:04d}-{1:02d}.npy"
NPY_3RD_FILE = "npydata-{0:04d}-{1:02d}-{2:02d}.npy"
L03B_2ND_FILE = "npyl03-b-{0:04d}-{1:02d}.npy"
L03B_3RD_FILE = "npyl03-b-{0:04d}-{1:02d}-{2:02d}.npy"

if len(sys.argv) < 3:
    sys.exit("Usage: loadnpy4.py mesh1 mesh2 [mesh3]")

mc = minecraft.Minecraft()

mesh1 = int(sys.argv[1])
mesh2 = int(sys.argv[2])
mesh3 = None
if len(sys.argv) >= 4:
    mesh3 = int(sys.argv[3])

if mesh3 is None:
    sys.exit("No implement !")

h_array = np.load(nd.NPY_DIR+"/"+NPY_3RD_FILE.format(mesh1, mesh2, mesh3))
l_array = np.load(nd.NPY_DIR+"/"+L03B_3RD_FILE.format(mesh1, mesh2, mesh3))

for z in range(150):
    for x in range(225):
        hval = h_array[z][x]
        lval = l_array[z][x]
        hval *= 0.2
        mc.setBlocks(x, hval-1, z, x, hval-7, z, block.DIRT)
        if lval == 100:  # 田
            mc.setBlock(x, hval, z, block.FARMLAND)
            mc.setBlock(x, hval+1, z, block.WHEAT)
        elif lval == 200: # その他農地
            mc.setBlock(x, hval, z, block.HARDENED_CLAY_STAINED_ORANGE)
        elif lval == 500: # 森林
            mc.setBlock(x, hval, z, block.GRASS)
            r = random.random()
            if r < 0.02:
                wood.wood(mc, x, hval+1, z)
            elif r < 0.1:
                mc.setBlock(x, hval+1, z, block.FERN)
        elif lval == 600: # 荒地
            mc.setBlock(x, hval, z, block.HARDENED_CLAY_STAINED_YELLOW)
        elif lval == 700: # 建物
            mc.setBlock(x, hval, z, block.HARDENED_CLAY_STAINED_RED)
        elif lval == 901: # 道路
            mc.setBlock(x, hval, z, block.HARDENED_CLAY_STAINED_LIGHT_GRAY)
        elif lval == 902: # 鉄道
            mc.setBlock(x, hval, z, block.HARDENED_CLAY_STAINED_GRAY)
        elif lval == 1000: # その他用地
            mc.setBlock(x, hval, z, block.HARDENED_CLAY_STAINED_PINK)
        elif lval == 1100: # 河川及び湖沼
            mc.setBlock(x, hval, z, block.HARDENED_CLAY_STAINED_BLUE)
        elif lval == 1400: # 海浜
            mc.setBlock(x, hval, z, block.HARDENED_CLAY_STAINED_CYAN)
        elif lval == 1500: # 海水域
            mc.setBlock(x, hval, z, block.HARDENED_CLAY_STAINED_PURPLE)
        elif lval == 1600: # ゴルフ場
            mc.setBlock(x, hval, z, block.HARDENED_CLAY_STAINED_LIME)
        else:
            mc.setBlock(x, hval, z, block.DIRT)
    sleep(0.02)

2017-10-10

標高データに土地利用細分メッシュを重ねる（１）

ジオクラ

標高データで地形を作りながら、一番上のブロックを土地利用細分メッシュで色分けしてみた。

f:id:pycra:20171010004232p:plain

いざ、塗ってみるとちょっと大雑把。

/python loadnpy4.py 5340 22 0

スクリプト

loadnpy4.pyとして保存。

# coding: utf-8

import mcpi.minecraft as minecraft
import mcpi.block as block
import numpy as np
import pycra.npydata as nd
import sys
from time import sleep

NPY_2ND_FILE = "npydata-{0:04d}-{1:02d}.npy"
NPY_3RD_FILE = "npydata-{0:04d}-{1:02d}-{2:02d}.npy"
L03B_2ND_FILE = "npyl03-b-{0:04d}-{1:02d}.npy"
L03B_3RD_FILE = "npyl03-b-{0:04d}-{1:02d}-{2:02d}.npy"

if len(sys.argv) < 3:
    sys.exit("Usage: loadnpy4.py mesh1 mesh2 [mesh3]")

mc = minecraft.Minecraft()

mesh1 = int(sys.argv[1])
mesh2 = int(sys.argv[2])
mesh3 = None
if len(sys.argv) >= 4:
    mesh3 = int(sys.argv[3])

if mesh3 is None:
    sys.exit("No implement !")

h_array = np.load(nd.NPY_DIR+"/"+NPY_3RD_FILE.format(mesh1, mesh2, mesh3))
l_array = np.load(nd.NPY_DIR+"/"+L03B_3RD_FILE.format(mesh1, mesh2, mesh3))

for z in range(150):
    for x in range(225):
        hval = h_array[z][x]
        lval = l_array[z][x]
        hval *= 0.2
        mc.setBlocks(x, hval-1, z, x, hval-7, z, block.DIRT)
        if lval == 100:  # 田
            mc.setBlock(x, hval, z, block.HARDENED_CLAY_STAINED_BROWN)
        elif lval == 200: # その他農地
            mc.setBlock(x, hval, z, block.HARDENED_CLAY_STAINED_ORANGE)
        elif lval == 500: # 森林
            mc.setBlock(x, hval, z, block.HARDENED_CLAY_STAINED_GREEN)
        elif lval == 600: # 荒地
            mc.setBlock(x, hval, z, block.HARDENED_CLAY_STAINED_YELLOW)
        elif lval == 700: # 建物
            mc.setBlock(x, hval, z, block.HARDENED_CLAY_STAINED_RED)
        elif lval == 901: # 道路
            mc.setBlock(x, hval, z, block.HARDENED_CLAY_STAINED_LIGHT_GRAY)
        elif lval == 902: # 鉄道
            mc.setBlock(x, hval, z, block.HARDENED_CLAY_STAINED_GRAY)
        elif lval == 1000: # その他用地
            mc.setBlock(x, hval, z, block.HARDENED_CLAY_STAINED_PINK)
        elif lval == 1100: # 河川及び湖沼
            mc.setBlock(x, hval, z, block.HARDENED_CLAY_STAINED_BLUE)
        elif lval == 1400: # 海浜
            mc.setBlock(x, hval, z, block.HARDENED_CLAY_STAINED_CYAN)
        elif lval == 1500: # 海水域
            mc.setBlock(x, hval, z, block.HARDENED_CLAY_STAINED_PURPLE)
        elif lval == 1600: # ゴルフ場
            mc.setBlock(x, hval, z, block.HARDENED_CLAY_STAINED_LIME)
        else:
            mc.setBlock(x, hval, z, block.DIRT)
    sleep(0.02)

2017-10-09

土地利用細分メッシュから切り出す

ジオクラライブラリ

土地利用細分メッシュから指定の範囲を切り出すスクリプト(detach.py)を書いた。

# coding: utf-8

import numpy as np
import matplotlib.pyplot as plot
import pycra.npydata as nd
import sys

L03B_FORMAT = "npyl03-b-{0:04d}.npy"
L03B_MESH2_FORMAT = "npyl03-b-{0:04d}-{1:02d}.npy"
L03B_MESH3_FORMAT = "npyl03-b-{0:04d}-{1:02d}-{2:02d}.npy"

if len(sys.argv) < 3:
    sys.exit("Usage: detach.py mesh1 mesh2 [mesh3]")

mesh1 = int(sys.argv[1])
mesh2 = int(sys.argv[2])

mesh3 = None
if len(sys.argv) >= 4:
    mesh3 = int(sys.argv[3])

array = np.load(nd.NPY_DIR+"/"+L03B_FORMAT.format(mesh1))

min_y = 0
min_x = 0
(max_y, max_x) = array.shape

max_y = max_y - (mesh2 // 10) * 100
min_y = max_y - 100
min_x = (mesh2 % 10) * 100
max_x = min_x + 100
print("min_y=", min_y)
print("max_y=", max_y)
print("min_x=", min_x)
print("max_x=", max_x)

if mesh3 is None:
    out_array = array[min_y:max_y, min_x:max_x]
    array15 = np.repeat(out_array, 15, axis=0)
    array15x23 = np.repeat(array15, 23, axis=1)
    del_collect = []
    for i in range(50):
        del_collect.append((i+1)*46-1)
    out_array = np.delete(array15x23, del_collect, 1)
    np.save(nd.NPY_DIR+"/"+L03B_MESH2_FORMAT.format(mesh1, mesh2), out_array)
else:
    max_yy = max_y - (mesh3 // 10) * 10
    min_yy = max_yy - 10
    min_xx = min_x + (mesh3 % 10) * 10
    max_xx = min_xx + 10
    print("min_yy=", min_yy)
    print("max_yy=", max_yy)
    print("min_xx=", min_xx)
    print("max_xx=", max_xx)
    out_array = array[min_yy:max_yy, min_xx:max_xx]
    print(out_array)
    array15 = np.repeat(out_array, 15, axis=0)
    array15x23 = np.repeat(array15, 23, axis=1)
    del_collect = []
    for i in range(5):
        del_collect.append((i+1)*46-1)
    out_array = np.delete(array15x23, del_collect, 1)
    np.save(nd.NPY_DIR+"/"+L03B_MESH3_FORMAT.format(mesh1, mesh2, mesh3), out_array)

print(out_array.shape)

plot.imshow(out_array, interpolation='nearest')
plot.colorbar()
plot.show()

実行結果

python detach.py 5340 22

f:id:pycra:20171009230312p:plain

python detach.py 5340 22 0

f:id:pycra:20171009230332p:plain

Pythonメモ

NumPyのarrayから行または列を指定して削除する。

keisanbutsuriya.hateblo.jp

repeatで23倍してからdeleteで間引いた。

2017-10-09

土地利用細分メッシュと基盤地図情報のメッシュの対応

国土数値情報基盤地図情報

土地利用細分メッシュと基盤地図情報の対応範囲を確認した。

土地利用細分メッシュは若干縦に伸びた感じ。

f:id:pycra:20171008122739p:plain

f:id:pycra:20171009110306p:plain

土地利用の方は3次メッシュとか書かれているけど、基盤地図情報の1次メッシュの範囲。

それぞれのメッシュの対応をまとめてみた。

f:id:pycra:20171012055745p:plain

土地利用細分メッシュの10x10が基盤地図情報の3次メッシュの150x225に対応している。

縦の10→150はぴったり15倍だからいいけど、横の10→225が問題。まずは22倍するのと23倍するのを交互にやって様子を見る。

2017-10-09

土地利用細分メッシュ再考

ジオクラ

土地利用細分メッシュが100mメッシュだからと言って、まるっきり使えないというのはもったいないので、「各種別毎に表示してみたらどうだろうか？」と種別コードを指定して表示するスクリプトを書いて実行してみた。

実行結果

個別に出力すると意外と使えそうなデータもありそう。

例えば、田、その他農用地とか、森林とかはどうせ広いからこの程度のメッシュでも使えるような気がする。

また、荒地、海浜とかは地表に石とか砂が露出したイメージで草、土ブロック一辺倒から変化をつけられる。

道路、鉄道は微妙。ちゃんとつながってないのが痛い。

種別100（田）

f:id:pycra:20171009053313p:plain

種別200（その他農用地）

f:id:pycra:20171009053327p:plain

種別500（森林）

f:id:pycra:20171009053342p:plain

種別600（荒地）

f:id:pycra:20171009053355p:plain

種別700（建物用地）

f:id:pycra:20171009053407p:plain

種別901（道路）

f:id:pycra:20171009053419p:plain

種別902（鉄道）

f:id:pycra:20171009053436p:plain

種別1000（その他用地）

運動競技場、空港、競馬場、野球場、学校、港湾地区、人工造成地の空地 f:id:pycra:20171009053449p:plain

種別1100（河川地及び湖沼）

f:id:pycra:20171009053501p:plain

種別1400（海浜）

f:id:pycra:20171009053514p:plain

種別1500（海水域）

f:id:pycra:20171009053524p:plain

種別1600（ゴルフ場）

f:id:pycra:20171009053535p:plain

スクリプト

# coding: utf-8

import xml.etree.ElementTree as ET
import numpy as np
import matplotlib.pyplot as plot
import sys

JPGIS_FILE = "JPGIS/5340/L03-b-14_5340-jgd_GML/L03-b-14_5340.xml"

if len(sys.argv) < 2:
    sys.exit("Usage: jgsb2npy.py typeNo")

typeNo = int(sys.argv[1])

tree = ET.parse(JPGIS_FILE)
root = tree.getroot()

for e in root.getiterator():
    print(e.tag)

tl = root.find('./{http://nlftp.mlit.go.jp/ksj/schemas/ksj-app}LanduseSubdivisionMesh/{http://nlftp.mlit.go.jp/ksj/schemas/ksj-app}coverage/{http://www.opengis.net/gml/3.2}rangeSet/{http://www.opengis.net/gml/3.2}DataBlock/{http://www.opengis.net/gml/3.2}tupleList')
if tl is None:
    raise Exception("{http://www.opengis.net/gml/3.2}tupleList is not found")

lines = tl.text.split("\n")
array = np.empty((800,800), dtype=int)
i = 0
for l in lines:
    d = l.split(" ")
    d.pop(0)
    d = np.array(d, dtype=int)
    if len(d) == 800:
        array[i] = d
        i += 1

array = 1 * (array == typeNo)

plot.imshow(array, interpolation='nearest')
plot.colorbar()
plot.show()

2017-10-08

土地利用細分メッシュ

ジオクラ

国土数値情報ダウンロードサービス

続いて、土地利用細分メッシュもマップとして表示してみた。

f:id:pycra:20171008122739p:plain

細分とはいえ、100mメッシュなのでこれも厳しい。

# coding: utf-8

import xml.etree.ElementTree as ET
import numpy as np
import matplotlib.pyplot as plot
import sys


JPGIS_FILE = "JPGIS/5340/L03-b-14_5340-jgd_GML/L03-b-14_5340.xml"

tree = ET.parse(JPGIS_FILE)
root = tree.getroot()

for e in root.getiterator():
    print(e.tag)

tl = root.find('./{http://nlftp.mlit.go.jp/ksj/schemas/ksj-app}LanduseSubdivisionMesh/{http://nlftp.mlit.go.jp/ksj/schemas/ksj-app}coverage/{http://www.opengis.net/gml/3.2}rangeSet/{http://www.opengis.net/gml/3.2}DataBlock/{http://www.opengis.net/gml/3.2}tupleList')
if tl is None:
    raise Exception("{http://www.opengis.net/gml/3.2}tupleList is not found")

lines = tl.text.split("\n")
print("lines=", len(lines))
array = np.empty((800,800), dtype=int)
i = 0
for l in lines:
    print("i=",i)
    d = l.split(" ")
    d.pop(0)
    d = np.array(d, dtype=int)
    print("len(d)=",len(d))
    if len(d) == 800:
        print(d)
        print(len(array[i]))
        array[i] = d
        i += 1

plot.imshow(array, interpolation='nearest')
plot.colorbar()
plot.show()

パイクラおじさんの日記

MinecraftでPythonを勉強するおじさんの日記です。

標高データに土地利用細分メッシュを重ねる（３）

スクリプト

標高データに土地利用細分メッシュを重ねる（２）

スクリプト

標高データに土地利用細分メッシュを重ねる（１）

スクリプト

土地利用細分メッシュから切り出す

実行結果

Pythonメモ

土地利用細分メッシュと基盤地図情報のメッシュの対応

土地利用細分メッシュ再考

実行結果

種別100（田）

種別200（その他農用地）

種別500（森林）

種別600（荒地）

種別700（建物用地）

種別901（道路）

種別902（鉄道）

種別1000（その他用地）

種別1100（河川地及び湖沼）

種別1400（海浜）

種別1500（海水域）

種別1600（ゴルフ場）

スクリプト

土地利用細分メッシュ