データ数が少ないエリア(1)
534023でデータ数が少ないエリアがあったので、それがどこなのかを調べたところ最初に見つかったのが04だった。
xmlplot.pyで処理してみると、33645行あって150 x 225 = 33750行あるべきなのに足りないのは…
と、スタート位置<gml:startPoint>を確認すると0 0で正しい。
なのでxmlplot.pyを修正して、最初から配列を33750行分割り当てて、頭から読み込ませて残ったらそのままにするようにした。
一番下の1行の右側に0mの青緑色の線が残っている。 ここは、海水面だと思うので、最初に全体を-9で埋めておけば良さそう。
配列を-9で埋めるのは.fill(値)でできる。
その修正をplotall.pyに入れて実行したのが、これ。
ちゃんと海岸線が表示された。多分、海岸線が左上にある場合にスタート位置<gml:startPoint>を確認する必要があるんだと思う。
5340-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()