パイクラおじさんの日記

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

データ数が少ないエリア(1)

534023でデータ数が少ないエリアがあったので、それがどこなのかを調べたところ最初に見つかったのが04だった。

xmlplot.pyで処理してみると、33645行あって150 x 225 = 33750行あるべきなのに足りないのは…

と、スタート位置<gml:startPoint>を確認すると0 0で正しい。

なのでxmlplot.pyを修正して、最初から配列を33750行分割り当てて、頭から読み込ませて残ったらそのままにするようにした。

f:id:pycra:20170918030359p:plain
Fig. 23-04

一番下の1行の右側に0mの青緑色の線が残っている。 ここは、海水面だと思うので、最初に全体を-9で埋めておけば良さそう。

配列を-9で埋めるのは.fill(値)でできる。

f:id:pycra:20170918031006p:plain
Fig. 23-04(改)

その修正をplotall.pyに入れて実行したのが、これ。

f:id:pycra:20170918031406p:plain
Fig. 23(改)

ちゃんと海岸線が表示された。多分、海岸線が左上にある場合にスタート位置<gml:startPoint>を確認する必要があるんだと思う。

5340-13部分もちゃんと海岸線が表示されるようになった。

f:id:pycra:20170918032141p:plain
Fig. 13(改)

現在のplotall.py。

# coding: utf-8

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

GEO_DIR = "FG-GML-{0:04d}-{1:02d}-DEM5A"
GEO_XML = "FG-GML-{0:04d}-{1:02d}-{2:02d}-DEM5A-20161001.xml"

def xml2array(mesh1, mesh2, mesh3):
    dir = GEO_DIR.format(mesh1, mesh2)
    fname = GEO_XML.format(mesh1, mesh2, mesh3)

    tree = ET.parse(dir+"/"+fname)
    root = tree.getroot()

    tl = root.find('./{http://fgd.gsi.go.jp/spec/2008/FGD_GMLSchema}DEM/{http://fgd.gsi.go.jp/spec/2008/FGD_GMLSchema}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()
    array = np.zeros(33750)
    array.fill(-9)
    i = 0
    for l in lines:
        (t, h) = l.split(",")
        hval = float(h)
        if hval == -9999:
            array[i] = -9
        else:
            array[i] = hval
        i += 1

    return array.reshape((150, 225))


# ここからスタート

if len(sys.argv) < 3:
    print("plotall.py mesh1 mesh2")
    exit(-1)

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

yarray = None
for y in range(10):
    xarray = None
    for x in range(10):
        try:
            mesh3 = y * 10 + x
            array = xml2array(mesh1, mesh2, mesh3)
        except Exception as err:
            print(err)
            array = np.zeros((150, 225))
            array.fill(-10)
        finally:
            if xarray is None:
                xarray = array
            else:
                xarray = np.concatenate((xarray, array), axis=1)
        
    if yarray is None:
        yarray = xarray
    else:
        yarray = np.concatenate((xarray, yarray), axis=0)

print(yarray.shape)

plot.imshow(yarray)
plot.colorbar()
plot.show()