立体地形図を手作りしてみようかな、っと(11.725)

(12)の続きではなくて、(11.5)の進化、のはなし。デクリメントしてどーする。

(11.5)で「メッシュ番号との対応がわかるとなお良いのだが」なんて書いたけど、よく考えれば、こんなん難しいわきゃねーな、と思い直し。

これ(basemap のソースツリーに含まれる examples)、なんだけど、プレインテキストで読みにくいので、全引用で pygmentize しとく:

customticks.py
 1 from __future__ import unicode_literals
 2 from mpl_toolkits.basemap import Basemap
 3 import numpy as np
 4 import matplotlib.pyplot as plt
 5 from matplotlib.ticker import FuncFormatter
 6 
 7 # example showing how to create custom tick labels for a cylindrical
 8 # projection.
 9 
10 def lat2str(deg):
11     min = 60 * (deg - np.floor(deg))
12     deg = np.floor(deg)
13     dir = 'N'
14     if deg < 0:
15         if min != 0.0:
16             deg += 1.0
17             min -= 60.0
18         dir = 'S'
19     return ("%d\N{DEGREE SIGN} %g' %s") % (np.abs(deg),np.abs(min),dir)
20 
21 def lon2str(deg):
22     min = 60 * (deg - np.floor(deg))
23     deg = np.floor(deg)
24     dir = 'E'
25     if deg < 0:
26         if min != 0.0:
27             deg += 1.0
28             min -= 60.0
29         dir = 'W'
30     return ("%d\N{DEGREE SIGN} %g' %s") % (np.abs(deg),np.abs(min),dir)
31 
32 # (1) use matplotlib custom tick formatter
33 #     instead of Basemap labelling methods.
34 
35 # create figure.
36 fig=plt.figure()
37 # create Basemap instance (regular lat/lon projection).
38 # suppress_ticks=False allows custom axes ticks to be used
39 # Ticks are suppressed by default, so Basemap methods
40 # drawparallels and drawmeridians used to draw labelled lat/lon grid.
41 m = Basemap(llcrnrlon=-156.5,llcrnrlat=18.75,urcrnrlon=-154.5,urcrnrlat=20.5,
42             resolution='h',projection='cyl',suppress_ticks=False)
43 # draw coastlines, fill land and lake areas.
44 m.drawcoastlines()
45 m.fillcontinents(color='coral',lake_color='aqua')
46 # background color will be used for oceans.
47 m.drawmapboundary(fill_color='aqua')
48 # get axes instance.
49 ax = plt.gca()
50 # add custom ticks.
51 # This only works for projection='cyl'.
52 def xformat(x, pos=None): return lon2str(x)
53 xformatter = FuncFormatter(xformat)
54 ax.xaxis.set_major_formatter(xformatter)
55 def yformat(y, pos=None): return lat2str(y)
56 yformatter = FuncFormatter(yformat)
57 ax.yaxis.set_major_formatter(yformatter)
58 ax.fmt_xdata = lambda x: lon2str(x)
59 ax.fmt_ydata = lambda y: lat2str(y)
60 ax.grid()
61 ax.set_title('Hawaii')
62 
63 # (2) use Basemap labelling methods, but pass a
64 #     custom formatting function with the 'fmt' keyword.
65 
66 # create figure.
67 fig = plt.figure()
68 # create Basemap instance.
69 m = Basemap(llcrnrlon=-156.5,llcrnrlat=18.75,urcrnrlon=-154.5,urcrnrlat=20.5,
70             resolution='h',projection='cyl')
71 # draw coastlines, fill land and lake areas.
72 m.drawcoastlines()
73 m.fillcontinents(color='coral',lake_color='aqua')
74 # background color will be used for oceans.
75 m.drawmapboundary(fill_color='aqua')
76 # label meridians and parallels, passing string formatting function
77 # with 'fmt' keyword.
78 m.drawparallels(np.linspace(18,21,7),labels=[1,0,0,0],fmt=lat2str,dashes=[2,2])
79 m.drawmeridians(np.linspace(-157,-154,7),labels=[0,0,0,1],fmt=lon2str,dashes=[2,2])
80 plt.title('Hawaii')
81 
82 plt.show()

つまり FuncFormatter を set_major_formatter する割と普通の方法と、drawparallels、drawmeridians に fmt を渡す方法の2つで可能だ、と言っている。

すわ、真似ぶべし:

done_vis.py
 1 # -*- coding: utf-8 -*-
 2 import os
 3 import re
 4 import numpy as np
 5 import np_utils
 6 
 7 
 8 # ##################################################
 9 #
10 import datafile_tree
11 
12 
13 _all = np.zeros((135 * 8 * 10, 80 * 8 * 10))
14 for fn, dty, m1, m2, m3 in datafile_tree.walk("__converted"):
15     y = m1[0] * 8 * 10 + m2[0] * 10 + m3[0]
16     x = m1[1] * 8 * 10 + m2[1] * 10 + m3[1]
17     print(dty, m1, m2, m3, (y, x))
18     _all[y, x] += 2. if dty == "DEM05" else 1.
19 
20 bd = np_utils.get_bound(_all, f=lambda a: a == 0.)
21 _x2lon = lambda x: x / (8 * 10.) + 100.
22 _y2lat = lambda y: y / (8 * 10.) * (2 / 3.)
23 elevs = _all[bd]
24 X = np.linspace(
25     _x2lon(bd[1].start),
26     _x2lon(bd[1].stop - 1),
27     elevs.shape[1])
28 Y = np.linspace(
29     _y2lat(bd[0].start),
30     _y2lat(bd[0].stop - 1),
31     elevs.shape[0])
32 
33 
34 # ##################################################
35 #
36 from mpl_toolkits.basemap import Basemap
37 import matplotlib.pyplot as plt
38 import matplotlib.cm as cm
39 
40 fig, ax = plt.subplots()
41 fig.set_size_inches(16.53 * 2, 11.69 * 2)
42 m = Basemap(
43     llcrnrlon=X[0] - (X[10] - X[0]) / 2.,
44     llcrnrlat=Y[0] - (Y[10] - Y[0]) / 2.,
45     urcrnrlon=X[-1] + (X[2 * 10] - X[0]),
46     urcrnrlat=Y[-1],
47     resolution='i',
48     lat_0=(Y[0] + Y[-1]) / 2.,
49     lon_0=(X[0] + X[-1]) / 2.,
50     projection='cass')
51 
52 lons, lats = np.meshgrid(X, Y)
53 x, y = m(lons, lats)
54 
55 m.drawcoastlines(color='0.9', linewidth=1)
56 #m.fillcontinents(alpha=0.3)
57 
58 def lon2str(deg):
59     return "%.1f\n[%d]" % (deg, deg - 100)
60 
61 def lat2str(deg):
62     return "%.3f\n[%d]" % (deg, deg / (2. / 3))
63 
64 circles = np.arange(10, 90, (2 / 3.))
65 m.drawparallels(circles, labels=[1, 0, 0, 0],
66                 fmt=lat2str)
67 meridians = np.arange(-180, 180, 1.)
68 m.drawmeridians(meridians, labels=[0, 0, 0, 1],
69                 fmt=lon2str)
70 #
71 m.pcolor(x, y, elevs, cmap=cm.summer)
72 plt.savefig("motterumotteru.png", bbox_inches="tight")
73 #plt.show()

基盤地図情報数値標高モデルのメッシュ番号が、緯度方向は「北半球を135分割(1次メッシュ)・さらに8分割(2次メッシュ)・さらに10分割(3次メッシュ)」、経度方向は「経度-100(一次メッシュ)・さらに8分割(2次メッシュ)・さらに10分割(3次メッシュ)」ね、念のためおさらい。ので parallels は 2/3 ずつのステップにすると、ちょうど1次メッシュの区切りになる。

結果はこんなだよ:

2次メッシュまでわかりやすくしても良さげではあるが…、かえって読みづらいかね?