やってみるだけやってみたならメモだけはしときたい、という話。
本当に簡単ならこんな書き方はしねー、ってか。
基本的にワタシのスタンスは、cartopy ネタの(5)で言ってみせたように、原則「必要とならない限りはミニマルで」である。ゆえに「GIS の一部が必要だ」のだとしても、「データベースと GIS の連携」は、今の自分には必要ないので、考えてない。そもそもワタシの最初の「GIS 体験」は PostGIS から始まったので、無論これを導入した経験もあるし、それどころか「かなりディープなアプリケーションを作っていた」。けど、今はそうではないし、個人ユースでそこまで必要になる状況になってない、ここ何年も。
ゆえに、まぁいつもながらの「皮算用」というか。全然自分では必要としてないんだけれど、「必要になったとしたら?」の話。「GIS のうち、ミニマルではなく、重量級」の一翼である「GIS とデータベース」の中でも、「SQLite なら「ライト」なんじゃねーの?」と思ってやってみた、て話。やってみたら確かに「考えようによってはライトだった」。
これから書く「答えの一つ」は、たぶんほとんどのユーザをがっかりさせるものである。けれども、少なくともワタシにとっては「なんつーチョロい」ものであるし、あるいは「野良ビルド上等」派の人にとっては結構な朗報なんではないのかな、という、かなり評価が割れるであろう情報ね。それを予め覚悟しといてちょ。
「簡単だ」と言えるのは、箇条書きで書けばこれだけだからさ:
- spatialite の、vcpkg でのビルドは失敗する。
- けれどもその依存物のビルドはすべて成功する。
- ゆえに、ソースコードからビルドしてしまえばいい。
- (64bit版の場合は)makefile_mod64.vc を書き換えてビルド(
nmake -f makefile_mod64.vc)すればいいだけ。C:/OSGeo4W64 への依存を vcpkg のものに書き換えればよろしい。ただし「_i.lib」になってるものは「_i」を取り除く。 - epsg_inlined_prussian.c のエラーはご随意に。使わないでしょ、ワタシは全部コメントアウトした。
- ビルド出来た mod_spatialite.dll は Python 3.6+ に同梱されてる sqlite3.dll と互換性がないので、Python3x/DLLs のものを vcpkg のものに差し替える。
- ビルド出来た mod_spatialite.dll は PATH が通ってる場所に置く。
- さすれば sqlite3 モジュールはそのままで mod_spatialite をロード出来る。
5. と 6. が、大多数のユーザを怯ませるんだとは思うけれど、正直 vcpkg のおかげでどんだけ楽か、って思うのよね、あたしゃ。考え方次第だよさ、て思うの。
ただそうは言ってもまぁ…、こういう苦労を好んでするくらいならさ、「Unix 使えよ」てのが正論だし、たぶん Windows であっても cygwin、msys2 を前提にすれば遥かに簡単なのだと思うよ。それでもどうしても Windows native な CPython でなければならない、という人は、試してみれば?
とにかく、ここまで出来てれば、あとは「Windows 版 CPython」という但し書きは関係なくなり…:
InitSpatialMetadata せずに出来ること。
1 >>> import sqlite3
2 >>>
3 >>> conn = sqlite3.connect("example.db")
4 >>> conn.enable_load_extension(True)
5 >>> c = conn.cursor()
6 >>> c.execute("SELECT load_extension('mod_spatialite')")
7 <sqlite3.Cursor object at 0x0000023539B10960>
8 >>> c.execute("SELECT IsValid(ST_Point(10, 20))")
9 <sqlite3.Cursor object at 0x0000023539B10960>
10 >>> c.fetchone()
11 (1,)
12 >>> c.execute("SELECT AsWKT(ST_Point(10, 20))")
13 <sqlite3.Cursor object at 0x0000023539B10960>
14 >>> c.fetchone()
15 ('POINT(10 20)',)
16 >>> c.execute("SELECT AsWKT(GeomFromText('LINESTRING (1 2, 3 4, 4 5)'))")
17 <sqlite3.Cursor object at 0x0000023539B10960>
18 >>> c.fetchone()
19 ('LINESTRING(1 2,3 4,4 5)',)
20 >>>
21 >>> conn.close()
22 >>>
InitSpatialMetadata。とても時間がかかるよ。
1 >>> import sqlite3
2 >>>
3 >>> conn = sqlite3.connect("example.db")
4 >>> conn.enable_load_extension(True)
5 >>> c = conn.cursor()
6 >>>
7 >>> # Create and populate SpatiaLite metadata tables.
8 >>> c.execute("SELECT InitSpatialMetadata()")
9 <sqlite3.Cursor object at 0x0000023539B10960>
10 >>> conn.close()
11 >>>
InitSpatialMetadataしたセッションを閉じたあとのインスタンスにて。
1 >>> import sqlite3
2 >>>
3 >>> conn = sqlite3.connect("example.db")
4 >>> conn.enable_load_extension(True)
5 >>> c = conn.cursor()
6 >>> c.execute("SELECT load_extension('mod_spatialite')")
7 <sqlite3.Cursor object at 0x0000026B91801960>
8 >>>
9 >>> # srid, auth_name, auth_srid, ref_sys_name, proj4text, srtext
10 >>> c.execute("SELECT srid, auth_name, ref_sys_name FROM spatial_ref_sys limit 10")
11 <sqlite3.Cursor object at 0x0000026B91801960>
12 >>> for _ in range(10):
13 ... print(c.fetchone())
14 ...
15 (-1, 'NONE', 'Undefined - Cartesian')
16 (0, 'NONE', 'Undefined - Geographic Long/Lat')
17 (2000, 'epsg', 'Anguilla 1957 / British West Indies Grid')
18 (2001, 'epsg', 'Antigua 1943 / British West Indies Grid')
19 (2002, 'epsg', 'Dominica 1945 / British West Indies Grid')
20 (2003, 'epsg', 'Grenada 1953 / British West Indies Grid')
21 (2004, 'epsg', 'Montserrat 1958 / British West Indies Grid')
22 (2005, 'epsg', 'St. Kitts 1955 / British West Indies Grid')
23 (2006, 'epsg', 'St. Lucia 1955 / British West Indies Grid')
24 (2007, 'epsg', 'St. Vincent 45 / British West Indies Grid')
25 >>> conn.close()
26 >>>
テーブルも作ってみる。
1 >>> import sqlite3
2 >>>
3 >>> conn = sqlite3.connect("example.db")
4 >>> conn.enable_load_extension(True)
5 >>> c = conn.cursor()
6 >>> c.execute("SELECT load_extension('mod_spatialite')")
7 <sqlite3.Cursor object at 0x0000016401D91960>
8 >>>
9 >>> # Creating a Geometry-type at the same time the corresponding
10 >>> # table is created isn't allowed. You always must first create
11 >>> # the table, then adding the Geometry-column in a second time
12 >>> # and as a separate step.
13 >>> c.execute("""
14 ... CREATE TABLE t_my_geoms (
15 ... id INTEGER NOT NULL
16 ... PRIMARY KEY AUTOINCREMENT,
17 ... name TEXT NOT NULL);""")
18 <sqlite3.Cursor object at 0x0000016401D91960>
19 >>>
20 >>> # 4326, 'POINT', 'XY' で「SRID=4326 の POINT 型、Dimension model が [X, Y]」
21 >>> # の意味。Dimension model では、結構良くつかわれる [X, Y, measur value] (XYM)
22 >>> # も選べる。Dimension model のあとに、Not Null かどうかをコントロール出来る
23 >>> # 値を渡せる。1 を渡せば Not Null となる。(省略は 0 なので、デフォルトで
24 >>> # Null 許容になる。)
25 >>> c.execute("""
26 ... SELECT AddGeometryColumn('t_my_geoms', 'pos', 4326, 'POINT', 'XY');
27 ... """)
28 <sqlite3.Cursor object at 0x0000016401D91960>
29 >>>
30 >>> # Take into account that you may want to use spatial indexes.
31 >>> c.execute("SELECT CreateSpatialIndex('t_my_geoms', 'pos')")
32 <sqlite3.Cursor object at 0x0000016401D91960>
33 >>> conn.close()
34 >>>
挿入して、検索して。
1 >>> import sqlite3
2 >>>
3 >>> conn = sqlite3.connect("example.db")
4 >>> conn.enable_load_extension(True)
5 >>> c = conn.cursor()
6 >>> c.execute("SELECT load_extension('mod_spatialite')")
7 <sqlite3.Cursor object at 0x00000269904A1960>
8 >>>
9 >>> # Now, let's insert records.
10 >>> # SetSRID(ST_Point(...), 4326) は MakePoint(..., 4326) でも可。
11 >>> c.execute("""
12 ... INSERT INTO t_my_geoms (name, pos) VALUES
13 ... ('01', SetSRID(ST_Point(141.35437400, 43.06197200), 4326)),
14 ... ('13', SetSRID(ST_Point(139.69171700, 35.68956800), 4326))
15 ... """)
16 <sqlite3.Cursor object at 0x00000269904A1960>
17 >>>
18 >>> conn.commit()
19 >>>
20 >>> c.execute("SELECT name, AsWKT(pos) FROM t_my_geoms")
21 <sqlite3.Cursor object at 0x00000269904A1960>
22 >>> c.fetchall()
23 [('01', 'POINT(141.354374 43.06197199999999)'), ('13', 'POINT(139.691717 35.689568)')]
24 >>>
25 >>> c.execute("SELECT name, AsKML(pos) FROM t_my_geoms limit 1")
26 <sqlite3.Cursor object at 0x000002BC72331960>
27 >>> c.fetchall()
28 [('01', '<Point><coordinates>141.354374,43.06197199999999</coordinates></Point>')]
29 >>>
30 >>> c.execute("SELECT name, AsWKT(pos, 6) FROM t_my_geoms ORDER BY pos DESC limit 1")
31 <sqlite3.Cursor object at 0x0000022EE5362960>
32 >>> c.fetchall()
33 [('13', 'POINT(139.691717 35.689568)')]
34 >>>
35 >>> conn.close()
36 >>>
なお、spatialite 側でせっかく geometry 型で定義していても、上で言った通り「sqlite3 モジュールはそのまま」なわけで、ゆえに SQL からの戻りの geometry 型を python 側での「geometry を表現するナニカ」に読み替えるのは自分でやらなきゃいけないし、insert だの update で geometry 型を渡したい場合も同じ。読み替えはたとえば:
1 >>> import sqlite3
2 >>> import shapely.wkt
3 >>> import shapely.wkb
4 >>> import shapely.geometry as sgeom
5 >>>
6 >>> conn = sqlite3.connect("example.db")
7 >>> conn.enable_load_extension(True)
8 >>> c = conn.cursor()
9 >>> c.execute("SELECT load_extension('mod_spatialite')")
10 <sqlite3.Cursor object at 0x0000019DC77FB650>
11 >>>
12 >>> c.execute("SELECT name, AsWKT(pos, 6) FROM t_my_geoms")
13 <sqlite3.Cursor object at 0x0000019DC77FB650>
14 >>> for row in c.fetchall():
15 ... print(row[0], shapely.wkt.loads(row[1]))
16 ...
17 01 POINT (141.354374 43.061972)
18 13 POINT (139.691717 35.689568)
19 >>> c.execute("SELECT name, AsBinary(pos) FROM t_my_geoms")
20 <sqlite3.Cursor object at 0x0000019DC77FB650>
21 >>> for row in c.fetchall():
22 ... print(row[0], shapely.wkb.loads(row[1]))
23 ...
24 01 POINT (141.354374 43.061972)
25 13 POINT (139.691717 35.689568)
26 >>> c.execute("""
27 ... SELECT name, AsWKT(pos, 6) FROM t_my_geoms
28 ... WHERE ST_Within(pos, ST_PolygonFromText(?))
29 ... """, (sgeom.box(141, 43, 142, 44).wkt,))
30 <sqlite3.Cursor object at 0x0000019DC77FB650>
31 >>> for row in c.fetchall():
32 ... print(row[0], shapely.wkt.loads(row[1]))
33 ...
34 01 POINT (141.354374 43.061972)
35 >>> conn.close()
36 >>>
このような明示的な読み替え以外の手段については、後ほど。
念のためもう一度繰り返しておくけれど、基本的にはまずは「Windows ねばならぬ」とか「Windows native な CPython ねばならぬのだ」というのが「自分が置かれている状況や立場による制約」でなく単なる強迫観念なのであるならば、そんなもの、捨てておしまいっ、てことね。ワタシがこうして『出来るよ』とやってみせてるのは、そういうあなた向けではないの。どうしてもやむなき事情により、「Windows ネイティブの CPython である必要がある」人か、もしくは何らかの検証そのものが目的の人向け。
「検証」の一つとして「複数環境での振る舞いの違いそのものに興味がある」というのがある。開発のステータスによっては「Unix 版の特定の機能が壊れている」みたいなことって皆無ではなくて、そうした場合に「Windows 版ではどうか? もしも Windows 版が正しい振る舞いをしているのならば、それをベースにテストを記述出来る」とかね。まぁこれを必要とする状況というのは、考えれば結構ある。
でもそれってやっぱりかなり特殊な人だと思うわけよ。正直言って、多くの人、特に「仕事で PC を使う人々」には、いい加減「Winodows が最も簡単である」という思い込みを捨て去ってほしいと願っている。Unix のほうが何億倍も簡単なのだぞ? (もちろん状況によることは否定はしない。定年後の爺さん婆さんにも Unix が簡単だ、なんてことはきっとないよ、そこまでは言ってない。)
…みたいな、一応言っておかねばならぬことは言ったうえで。ことワタシに関しては、そうした事情はそれほどはなくて、強いて言えば、「多くの「Windows ねばならぬ事情の人々」」のために Windows ユーザであることを捨てない、というのが唯一の「制約」で、理由の本当のところは「使える Unix 実機環境は(このブログ運用してるレンタルサーバを除き)持ってない」「virtual box とかで Unix を動かすだとかがめんどい(「重い」のが一番の理由)」「msys2 や cygwin はデカ過ぎて、躊躇しがち」とかまぁそんなくだらない理由でもって「普段使いの公式 Windows native CPython で使えるのが一番便利」みたいな、まぁ消極的な理由だったりする。人に言っといてなんだよ、てことなんだけどさ。
まぁちょっと哲学とでも言うべき話だったかもしれんわけよ。んでもって、そうしたことなんかどうでもいい、という人でこの話をどう受け取ればいいのか悩むのならば、上述の導入手順がもたらす「管理・保守問題」だけに着目すればいいと思う。公式 CPython をインストールして置かれた「sqlite3.dll」を差し替えたりとかって、少なくともちゃんとメモしておかないと間違いなくそれをしたことを忘れてしまって、Python 本体をアップデートする際に上書きされてしまったりとか、そういうことが簡単に起こる。「非公式なことをする」ってのはそういうこと。
なお、「毎度 enable_load_extension して load_extension するのがダル過ぎる」に応えてくれる spatialite というややおバカめなライブラリがある。ほんとに少しのことしかやってくれないものだし、なので、全然なくて困るもんでもないけれど、結構楽は楽なので、まぁ使ってみてもいいんじゃないかな:
1 >>> import spatialite
2 >>>
3 >>> with spatialite.connect("example.db") as conn:
4 ... c = conn.cursor()
5 ... c.execute("SELECT name, AsWKT(pos) FROM t_my_geoms")
6 ... print(c.fetchall())
7 ...
8 <sqlite3.Cursor object at 0x000001B9A50D1CE0>
9 [('01', 'POINT(141.354374 43.06197199999999)'), ('13', 'POINT(139.691717 35.689568)')]
InitSpatialMetadata してくれるメソッドも提供されているので、たとえば:
1 >>> import spatialite
2 >>>
3 >>> with spatialite.connect("my_gis_memo.db") as conn:
4 ... conn.initialize_metadata()
5 ... c = conn.cursor()
6 ... c.execute("""
7 ... CREATE TABLE t_geoshape_city (
8 ... entry_id TEXT NOT NULL PRIMARY KEY,
9 ... body TEXT NOT NULL,
10 ... body_variants TEXT NOT NULL,
11 ... suffix TEXT NOT NULL,
12 ... prefname TEXT NOT NULL,
13 ... prefname_variants TEXT NOT NULL,
14 ... countyname TEXT NOT NULL,
15 ... countyname_variants TEXT NOT NULL,
16 ... ne_class TEXT NOT NULL,
17 ... hypernym TEXT NOT NULL,
18 ... address TEXT NOT NULL,
19 ... code TEXT NOT NULL,
20 ... valid_from DATE,
21 ... valid_to DATE,
22 ... source TEXT NOT NULL);""")
23 ... c.execute("""
24 ... SELECT AddGeometryColumn(
25 ... 't_geoshape_city', 'refpoint', 4326, 'POINT', 'XY', 1);""")
26 ... c.execute("SELECT CreateSpatialIndex('t_geoshape_city', 'refpoint')")
27 ...
28 <sqlite3.Cursor object at 0x000001EE659561F0>
29 <sqlite3.Cursor object at 0x000001EE659561F0>
30 <sqlite3.Cursor object at 0x000001EE659561F0>
InitSpatialMetadataの時間がかかることを考えると、一気にやろうとしない方がほんとは無難だがあくまで例として。
というわけで、geoshape-city の情報を丸ごと spatialite にブチこんじゃるぜ ver 0.0.0.1:
1 # -*- coding: utf-8 -*-
2 import io
3 import csv
4 import os
5 import itertools
6 from datetime import date
7 import shapely.geometry as sgeom
8 import spatialite
9
10
11 _DBNAME = "my_gis_memo.db"
12
13
14 def _initialize_db(conn):
15 conn.initialize_metadata()
16 c = conn.cursor()
17 c.execute("""
18 CREATE TABLE t_geoshape_city (
19 entry_id TEXT NOT NULL PRIMARY KEY,
20 body TEXT NOT NULL,
21 body_variants TEXT NOT NULL,
22 suffix TEXT NOT NULL,
23 prefname TEXT NOT NULL,
24 prefname_variants TEXT NOT NULL,
25 countyname TEXT NOT NULL,
26 countyname_variants TEXT NOT NULL,
27 ne_class TEXT NOT NULL,
28 hypernym TEXT NOT NULL,
29 address TEXT NOT NULL,
30 code TEXT NOT NULL,
31 valid_from DATE,
32 valid_to DATE,
33 source TEXT NOT NULL);""")
34 c.execute("""
35 SELECT AddGeometryColumn(
36 't_geoshape_city', 'refpoint', 4326, 'POINT', 'XY', 1);""")
37 c.execute("SELECT CreateSpatialIndex('t_geoshape_city', 'refpoint')")
38 conn.commit()
39
40
41 def _load_from_geoshape_city_csv():
42 # 『歴史的行政区域データセットβ版』（CODH作成）
43 # https://geonlp.ex.nii.ac.jp/dictionary/geoshape-city/
44 def _map(rec, fields):
45 import re
46 temp = []
47 for v in rec:
48 if re.match(r"\d+\.\d+", v):
49 temp.append(float(v))
50 elif re.match(r"\d+\-\d+\-\d+", v):
51 temp.append(date(*map(int, v.split("-"))))
52 else:
53 temp.append(v)
54 return {k: v for k, v in zip(fields, temp)}
55 reader = csv.reader(
56 io.open("geoshape-city.csv", encoding="utf-8"))
57 fields = next(reader)
58 return (_map(rec, fields) for rec in reader)
59
60
61 def _fill_data(conn):
62 fields = None
63 values = []
64 for rec in _load_from_geoshape_city_csv():
65 # 「代表点」
66 refpt = rec.pop("longitude"), rec.pop("latitude")
67 if not fields:
68 fields = list(sorted(rec.keys())) + ["refpoint", "refpoint_lat"]
69 values.append([rec[k] for k in fields[:-2]] + list(refpt))
70 #
71 c = conn.cursor()
72 chunk = 500
73 for i in range(0, len(values), chunk):
74 c.execute("""
75 INSERT INTO t_geoshape_city
76 ({})
77 VALUES
78 {}
79 """.format(
80 ", ".join(fields[:-1]),
81 ", ".join(["({})".format(
82 ", ".join(["?"] * len(row[:-2]) + ["MakePoint(?, ?, 4326)"]))
83 for row in values[i:i + chunk]])),
84 tuple(itertools.chain.from_iterable(values[i:i + chunk])))
85 conn.commit()
86
87
88 def _main(bb):
89 first = not os.path.exists(_DBNAME)
90 with spatialite.connect(_DBNAME) as conn:
91 if first:
92 _initialize_db(conn)
93 _fill_data(conn)
94 c = conn.cursor()
95 c.execute("""
96 SELECT address, code FROM t_geoshape_city
97 where ST_Within(refpoint, ST_PolyFromText(?))""", (sgeom.box(*bb).wkt,))
98 for found in c.fetchall():
99 print(found)
100
101
102 if __name__ == '__main__':
103 import json
104 import argparse
105 ap = argparse.ArgumentParser()
106 ap.add_argument("bb")
107 args = ap.parse_args()
108 _main(json.loads(args.bb))
1 [me@host: ~]$ ls -l geoshape-city.csv
2 -rw-r--r-- 3 hhsprings Administrators 4485002 Dec 3 2020 geoshape-city.csv
3 [me@host: ~]$ py -3 my_gis_memo '[123, 24, 124, 25]'
4 ('沖縄県八重山郡与那国町', '47382')
5 ('沖縄県八重山郡竹富村', '47381')
6 ('沖縄県八重山郡与那国村', '47382')
まったく苦労しなかったぜ、なんてのはウソである。「?」へのバインドの部分は頭痛くなるよね。こういう「ふたつでひとつ」を「ふたつ」として扱ったり「ひとつ」として扱ったりが場所によってまちまちになるのがめんどうでなぁ。(values の値としてはフラットに lon, lat と「ふたつ」でないと困り、「lon と lat はひとまとまり」と見る必要がある場所もある、てこと。)
「sqlite3.register_adapter」を使えば、「ふたつでひとつだったりふたつだったり」問題は少しだけ改善する:
1 # -*- coding: utf-8 -*-
2 import io
3 import csv
4 import os
5 import itertools
6 import sqlite3
7 from datetime import date
8 import shapely.geometry as sgeom
9 import spatialite
10
11
12 _DBNAME = "my_gis_memo.db"
13
14
15 def _initialize_db(conn):
16 conn.initialize_metadata()
17 c = conn.cursor()
18 c.execute("""
19 CREATE TABLE t_geoshape_city (
20 entry_id TEXT NOT NULL PRIMARY KEY,
21 body TEXT NOT NULL,
22 body_variants TEXT NOT NULL,
23 suffix TEXT NOT NULL,
24 prefname TEXT NOT NULL,
25 prefname_variants TEXT NOT NULL,
26 countyname TEXT NOT NULL,
27 countyname_variants TEXT NOT NULL,
28 ne_class TEXT NOT NULL,
29 hypernym TEXT NOT NULL,
30 address TEXT NOT NULL,
31 code TEXT NOT NULL,
32 valid_from DATE,
33 valid_to DATE,
34 source TEXT NOT NULL);""")
35 c.execute("""
36 SELECT AddGeometryColumn(
37 't_geoshape_city', 'refpoint', 4326, 'POINT', 'XY', 1);""")
38 c.execute("SELECT CreateSpatialIndex('t_geoshape_city', 'refpoint')")
39 conn.commit()
40
41
42 def _load_from_geoshape_city_csv():
43 # 『歴史的行政区域データセットβ版』（CODH作成）
44 # https://geonlp.ex.nii.ac.jp/dictionary/geoshape-city/
45 def _map(rec, fields):
46 import re
47 temp = []
48 for v in rec:
49 if re.match(r"\d+\.\d+", v):
50 temp.append(float(v))
51 elif re.match(r"\d+\-\d+\-\d+", v):
52 temp.append(date(*map(int, v.split("-"))))
53 else:
54 temp.append(v)
55 return {k: v for k, v in zip(fields, temp)}
56 reader = csv.reader(
57 io.open("geoshape-city.csv", encoding="utf-8"))
58 fields = next(reader)
59 return (_map(rec, fields) for rec in reader)
60
61
62 def _point_adapter(pt):
63 return pt.wkb
64 sqlite3.register_adapter(sgeom.Point, _point_adapter)
65
66
67 def _fill_data(conn):
68 fields = None
69 values = []
70 for rec in _load_from_geoshape_city_csv():
71 # 「代表点」
72 refpt = sgeom.Point(rec.pop("longitude"), rec.pop("latitude"))
73 if not fields:
74 fields = list(sorted(rec.keys())) + ["refpoint"]
75 values.append([rec[k] for k in fields[:-1]] + [refpt])
76 #
77 c = conn.cursor()
78 chunk = 500
79 for i in range(0, len(values), chunk):
80 c.execute("""
81 INSERT INTO t_geoshape_city
82 ({})
83 VALUES
84 {}
85 """.format(
86 ", ".join(fields),
87 ", ".join(["({})".format(
88 ", ".join(["?"] * len(row[:-1]) + ["GeomFromWKB(?, 4326)"]))
89 for row in values[i:i + chunk]])),
90 tuple(itertools.chain.from_iterable(values[i:i + chunk])))
91 conn.commit()
92
93
94 def _main(bb):
95 first = not os.path.exists(_DBNAME)
96 with spatialite.connect(_DBNAME) as conn:
97 if first:
98 _initialize_db(conn)
99 _fill_data(conn)
100 c = conn.cursor()
101 c.execute("""
102 SELECT address, code FROM t_geoshape_city
103 where ST_Within(refpoint, ST_PolyFromText(?))""", (sgeom.box(*bb).wkt,))
104 for found in c.fetchall():
105 print(found)
106
107
108 if __name__ == '__main__':
109 import json
110 import argparse
111 ap = argparse.ArgumentParser()
112 ap.add_argument("bb")
113 args = ap.parse_args()
114 _main(json.loads(args.bb))
shapely.geometry.Point という値が「一つの値である」という扱いを spatialite 側にまで伝播させる、的なことなのだけれど、説明しづらい…。「GeomFromWKB」の「?」一個として「Point#wkb」の値を渡せる、てこと。GeomFromWKB すらしない手段があればいいんだけど? WKB でもなく、EWKB ですらなくて弱っている。
「GeomFromWKB の必要」は「Geometry 型」というスペシャルなマッピングを経由していないバイナリの処し方、てこと。探したら、やっと見つけた。BLOB。39バイト目にバイトオーダ指示が来ていなくて38バイト目が「0x7c(“|”)」、ということなので、38バイト目から渡しちゃいけないんじゃないかと思ったら、これで読めるみたいなのよね:
1 c = conn.cursor()
2 # 「refpoint」が general BLOB にブチ込んでるだけなので、という返却をしでかす。
3 # (geometry specific な振る舞いは「頼まれるまではしない」というノリよね。)
4 c.execute("""
5 SELECT address, code, refpoint FROM t_geoshape_city
6 where ST_Within(refpoint, ST_PolyFromText(?))""", (sgeom.box(*bb).wkt,))
7 for found in c.fetchall():
8 blob = found[-1]
9 # 末尾が 0xfe (end mark) なので :-1 はいいとして、38: が腑に落ちない。
10 # けど、読めることは読める。
11 print(blob[38], "%x"%blob[38], blob[38:-1], shapely.wkb.loads(blob[38:-1]))
ともあれこの知識を用いると、sqlite3.register_converter を使って:
1 import shapely.wkb
2 import shapely.wkt
3
4 # ...
5
6 def _point_adapter(pt):
7 return pt.wkb
8 sqlite3.register_adapter(sgeom.Point, _point_adapter)
9 def _geomblob_converter(blob):
10 return shapely.wkb.loads(blob[38:-1])
11 sqlite3.register_converter("POINT", _geomblob_converter)
12
13 # ...
14 def _main(bb):
15 first = not os.path.exists(_DBNAME)
16 # sqlite3.PARSE_DECLTYPES 指定がないと register_converter で登録した
17 # ものは POINT 型にディスパッチされない。
18 with spatialite.connect(_DBNAME, detect_types=sqlite3.PARSE_DECLTYPES) as conn:
19 if first:
20 _initialize_db(conn)
21 _fill_data(conn)
22 c = conn.cursor()
23 c.execute("""
24 SELECT address, code, refpoint FROM t_geoshape_city
25 where ST_Within(refpoint, ST_PolyFromText(?))""", (sgeom.box(*bb).wkt,))
26 for found in c.fetchall():
27 print(found)
28
29 # ...
これで例えば SELECT * したとしても、ちゃんと POINT 型が shapely.geometry な型に読み替えられる。戻りはこれでいいが、行きの方な。この仕様に従ってエンコードするしかない。うーん、これ、「spatialite 固有」でないならいいんだけどなぁ? PostGIS も同じ、とかだったらハッピーなんだけれども、違う場合管理が面倒なんだよなぁ。どうしたものか…。とにかく、気持ちがあまり良くないコードだが、これで AsBinary だので明示的に読み替えなくてよくなる:
AsBinary とか ST_Point とかの明示的な読み替えを都度するのではなく…
1 # -*- coding: utf-8 -*-
2 import io
3 import csv
4 import os
5 import itertools
6 import struct
7 import sqlite3
8 from datetime import date
9 import shapely.geometry as sgeom
10 import shapely.wkb
11 import shapely.wkt
12 import spatialite
13
14
15 _DBNAME = "my_gis_memo.db"
16
17
18 def _initialize_db(conn):
19 conn.initialize_metadata()
20 c = conn.cursor()
21 c.execute("""
22 CREATE TABLE t_geoshape_city (
23 entry_id TEXT NOT NULL PRIMARY KEY,
24 body TEXT NOT NULL,
25 body_variants TEXT NOT NULL,
26 suffix TEXT NOT NULL,
27 prefname TEXT NOT NULL,
28 prefname_variants TEXT NOT NULL,
29 countyname TEXT NOT NULL,
30 countyname_variants TEXT NOT NULL,
31 ne_class TEXT NOT NULL,
32 hypernym TEXT NOT NULL,
33 address TEXT NOT NULL,
34 code TEXT NOT NULL,
35 valid_from DATE,
36 valid_to DATE,
37 source TEXT NOT NULL);""")
38 c.execute("""
39 SELECT AddGeometryColumn(
40 't_geoshape_city', 'refpoint', 4326, 'POINT', 'XY', 1);""")
41 c.execute("SELECT CreateSpatialIndex('t_geoshape_city', 'refpoint')")
42 conn.commit()
43
44
45 def _load_from_geoshape_city_csv():
46 # 『歴史的行政区域データセットβ版』（CODH作成）
47 # https://geonlp.ex.nii.ac.jp/dictionary/geoshape-city/
48 def _map(rec, fields):
49 import re
50 temp = []
51 for v in rec:
52 if re.match(r"\d+\.\d+", v):
53 temp.append(float(v))
54 elif re.match(r"\d+\-\d+\-\d+", v):
55 temp.append(date(*map(int, v.split("-"))))
56 else:
57 temp.append(v)
58 return {k: v for k, v in zip(fields, temp)}
59 reader = csv.reader(
60 io.open("geoshape-city.csv", encoding="utf-8"))
61 fields = next(reader)
62 return (_map(rec, fields) for rec in reader)
63
64
65 def _geom_adapter(g):
66 wkb = g.wkb
67 endian = wkb[0] # 0: big, 1: little
68 em = "<" if endian else ">"
69 bb = g.bounds
70 encoded = struct.pack(
71 em + "BBiddddB{}sB".format(len(wkb[1:])),
72 0x00,
73 endian,
74 4326, # SRID
75 bb[0], bb[1], bb[2], bb[3], # MBR
76 0x7c,
77 wkb[1:],
78 0xfe)
79 return encoded
80 sqlite3.register_adapter(sgeom.Point, _geom_adapter)
81 sqlite3.register_adapter(sgeom.Polygon, _geom_adapter)
82
83
84 def _geomblob_converter(blob):
85 return shapely.wkb.loads(blob[38:-1])
86 sqlite3.register_converter("POINT", _geomblob_converter)
87
88
89 def _fill_data(conn):
90 fields = None
91 values = []
92 for rec in _load_from_geoshape_city_csv():
93 # 「代表点」
94 refpt = sgeom.Point(rec.pop("longitude"), rec.pop("latitude"))
95 if not fields:
96 fields = list(sorted(rec.keys())) + ["refpoint"]
97 values.append([rec[k] for k in fields[:-1]] + [refpt])
98 #
99 c = conn.cursor()
100 chunk = 500
101 for i in range(0, len(values), chunk):
102 c.execute("""
103 INSERT INTO t_geoshape_city
104 ({})
105 VALUES
106 {}
107 """.format(
108 ", ".join(fields),
109 ", ".join(["({})".format(
110 ", ".join(["?"] * len(row)))
111 for row in values[i:i + chunk]])),
112 tuple(itertools.chain.from_iterable(values[i:i + chunk])))
113 conn.commit()
114
115
116 def _main(bb):
117 first = not os.path.exists(_DBNAME)
118 with spatialite.connect(_DBNAME, detect_types=sqlite3.PARSE_DECLTYPES) as conn:
119 if first:
120 _initialize_db(conn)
121 _fill_data(conn)
122 c = conn.cursor()
123 c.execute("""
124 SELECT address, code, refpoint FROM t_geoshape_city
125 where ST_Within(refpoint, ?)""", (sgeom.box(*bb),))
126 for found in c.fetchall():
127 print(found[:-1], found[-1].wkt)
128
129
130 if __name__ == '__main__':
131 import json
132 import argparse
133 ap = argparse.ArgumentParser()
134 ap.add_argument("bb")
135 args = ap.parse_args()
136 _main(json.loads(args.bb))
技術的なネタとしては以上で一段落ね。けどね、ここまで出来るとさ、当然欲が出てくるわけだ。駅データ.jpのデータも入ってると便利かしらん、と:
python 2.7 では動かないよ
1 # -*- coding: utf-8 -*-
2 import io
3 import csv
4 import os
5 import itertools
6 import struct
7 import sqlite3
8 from datetime import date
9 import shapely.geometry as sgeom
10 import shapely.wkb
11 import shapely.wkt
12 import spatialite
13
14
15 _DBNAME = "my_gis_memo.db"
16
17
18 #
19 #
20 #
21 def _geom_adapter(g):
22 wkb = g.wkb
23 endian = wkb[0] # 0: big, 1: little
24 em = "<" if endian else ">"
25 bb = g.bounds
26 encoded = struct.pack(
27 em + "BBiddddB{}sB".format(len(wkb[1:])),
28 0x00,
29 endian,
30 4326, # SRID
31 bb[0], bb[1], bb[2], bb[3], # MBR
32 0x7c,
33 wkb[1:],
34 0xfe)
35 return encoded
36 sqlite3.register_adapter(sgeom.Point, _geom_adapter)
37 sqlite3.register_adapter(sgeom.Polygon, _geom_adapter)
38
39
40 def _geomblob_converter(blob):
41 return shapely.wkb.loads(blob[38:-1])
42 sqlite3.register_converter("POINT", _geomblob_converter)
43 #
44 #
45 #
46
47
48 def _initialize_db(conn):
49 conn.initialize_metadata()
50 c = conn.cursor()
51
52 # -----------------------------------------------
53 # 『歴史的行政区域データセットβ版』（CODH作成）
54 c.execute("""
55 CREATE TABLE t_geoshape_city (
56 entry_id TEXT NOT NULL PRIMARY KEY,
57 body TEXT NOT NULL,
58 body_variants TEXT NOT NULL,
59 suffix TEXT NOT NULL,
60 prefname TEXT NOT NULL,
61 prefname_variants TEXT NOT NULL,
62 countyname TEXT NOT NULL,
63 countyname_variants TEXT NOT NULL,
64 ne_class TEXT NOT NULL,
65 hypernym TEXT NOT NULL,
66 address TEXT NOT NULL,
67 code TEXT NOT NULL,
68 valid_from DATE,
69 valid_to DATE,
70 source TEXT NOT NULL);""")
71 c.execute("""
72 SELECT AddGeometryColumn(
73 't_geoshape_city', 'refpoint', 4326, 'POINT', 'XY', 1);""")
74 c.execute("SELECT CreateSpatialIndex('t_geoshape_city', 'refpoint')")
75 #
76 # -----------------------------------------------
77
78 # -----------------------------------------------
79 # 『駅データ.jp』
80 #
81 # company
82 c.execute("""
83 CREATE TABLE t_ekidata_company (
84 company_cd INTEGER PRIMARY KEY NOT NULL,
85 rr_cd INTEGER NOT NULL,
86 company_name TEXT NOT NULL,
87 company_name_k TEXT NOT NULL,
88 company_name_h TEXT NOT NULL,
89 company_name_r TEXT NOT NULL,
90 company_url TEXT NOT NULL,
91 company_type INTEGER NOT NULL,
92 e_status INTEGER NOT NULL
93 )""")
94
95 # join
96 c.execute("""
97 CREATE TABLE t_ekidata_join (
98 line_cd INTEGER NOT NULL,
99 station_cd1 INTEGER NOT NULL,
100 station_cd2 INTEGER NOT NULL
101 )""")
102
103 # line
104 # lon, lat → refpoint
105 c.execute("""
106 CREATE TABLE t_ekidata_line (
107 line_cd INTEGER PRIMARY KEY NOT NULL,
108 company_cd INTEGER NOT NULL,
109 line_name TEXT NOT NULL,
110 line_name_k TEXT NOT NULL,
111 line_name_h TEXT NOT NULL,
112 e_status INTEGER NOT NULL
113 )""")
114 c.execute("""
115 SELECT AddGeometryColumn(
116 't_ekidata_line', 'refpoint', 4326, 'POINT', 'XY', 1)""")
117 c.execute("""
118 SELECT CreateSpatialIndex('t_ekidata_line', 'refpoint')""")
119
120 # station
121 # lon, lat → refpoint
122 c.execute("""
123 CREATE TABLE t_ekidata_station (
124 station_cd INTEGER PRIMARY KEY NOT NULL,
125 station_g_cd INTEGER NOT NULL,
126 station_name TEXT NOT NULL,
127 line_cd INTEGER NOT NULL,
128 pref_cd INTEGER NOT NULL,
129 post TEXT NOT NULL,
130 address TEXT NOT NULL,
131 open_ymd DATE,
132 close_ymd DATE,
133 e_status INTEGER NOT NULL
134 )""")
135 c.execute("""
136 SELECT AddGeometryColumn(
137 't_ekidata_station', 'refpoint', 4326, 'POINT', 'XY', 1)""")
138 c.execute("""
139 SELECT CreateSpatialIndex('t_ekidata_station', 'refpoint')""")
140 #
141 # ---------------------------------------------------
142 conn.commit()
143
144
145 def _load_from_geoshape_city_csv():
146 # 『歴史的行政区域データセットβ版』（CODH作成）
147 # https://geonlp.ex.nii.ac.jp/dictionary/geoshape-city/
148 def _map(rec, fields):
149 import re
150 temp = []
151 for v in rec:
152 if re.match(r"\d+\.\d+", v):
153 temp.append(float(v))
154 elif re.match(r"\d+\-\d+\-\d+", v):
155 temp.append(date(*map(int, v.split("-"))))
156 else:
157 temp.append(v)
158 return {k: v for k, v in zip(fields, temp)}
159 reader = csv.reader(
160 io.open("geoshape-city.csv", encoding="utf-8"))
161 fields = next(reader)
162 return (_map(rec, fields) for rec in reader)
163
164
165 # ----------------------------- 駅データ.jp の列定義
166 #
167 # company
168 coldefs_company = [
169 ("company_cd", {"ctor": int,}),
170 ("rr_cd", {"ctor": int,}),
171 ("company_name", {"ctor": str,}),
172 ("company_name_k", {"ctor": str,}),
173 ("company_name_h", {"ctor": str,}),
174 ("company_name_r", {"ctor": str,}),
175 ("company_url", {"ctor": str,}),
176 ("company_type", {"ctor": int,}),
177 ("e_status", {"ctor": int,}),
178 #("e_sort", {"ctor": int,}),
179 ]
180
181 # join
182 coldefs_join = [
183 ("line_cd", {"ctor": int,}),
184 ("station_cd1", {"ctor": int,}),
185 ("station_cd2", {"ctor": int,}),
186 ]
187
188 # line
189 coldefs_line = [
190 ("line_cd", {"ctor": int,}),
191 ("company_cd", {"ctor": int,}),
192 ("line_name", {"ctor": str,}),
193 ("line_name_k", {"ctor": str,}),
194 ("line_name_h", {"ctor": str,}),
195 #("line_color_c", {"ctor": int,}),
196 #("line_color_t", {"ctor": int,}),
197 #("line_type", {"ctor": int,}),
198 ("lon", {"ctor": float,}),
199 ("lat", {"ctor": float,}),
200 #("zoom", {"ctor": int,}),
201 ("e_status", {"ctor": int,}),
202 #("e_sort", {"ctor": int,}),
203 ]
204
205 # station
206 def _str2date(s):
207 # 0000-00-00 を無効値としているが、これは扱いやすいとはいえない。
208 # null としてしまうのが扱いやすいと思う。
209 tt = list(map(int, s.split("-")))
210 if sum(tt) > 0:
211 return date(*tt)
212
213
214 coldefs_station = [
215 ("station_cd", {"ctor": int,}),
216 ("station_g_cd", {"ctor": int,}),
217 ("station_name", {"ctor": str,}),
218 #("station_name_k", {"ctor": str,}),
219 #("station_name_r", {"ctor": str,}),
220 ("line_cd", {"ctor": int,}),
221 ("pref_cd", {"ctor": int,}),
222 ("post", {"ctor": str,}),
223 ("address", {"ctor": str,}),
224 ("lon", {"ctor": float,}),
225 ("lat", {"ctor": float,}),
226 ("open_ymd", {"ctor": _str2date,}),
227 ("close_ymd", {"ctor": _str2date,}),
228 ("e_status", {"ctor": int,}),
229 #("e_sort", {"ctor": int,}),
230 ]
231
232
233 def _load_from_ekidata_csv(csvfn, coldefs):
234 _DBKEYWORDS_TRANS = {
235 "add": "address",
236 }
237 with io.open(csvfn, encoding="utf-8") as fi:
238 reader = csv.reader(fi)
239 headers = next(reader)
240 headers = [_DBKEYWORDS_TRANS.get(v, v) for v in headers]
241 for vals in reader:
242 yield {k: ty["ctor"](vals[headers.index(k)])
243 for k, ty in coldefs}
244
245
246 def _fill_data(conn):
247 def _do_insert(cur, tblname, fields, records):
248 chunk = 500
249 for i in range(0, len(records), chunk):
250 cur.execute("""
251 INSERT INTO {}
252 ({})
253 VALUES
254 {}
255 """.format(
256 tblname,
257 ", ".join(fields),
258 ", ".join(["({})".format(
259 ", ".join(["?"] * len(row)))
260 for row in records[i:i + chunk]])),
261 tuple(itertools.chain.from_iterable(records[i:i + chunk])))
262 #
263 # 『歴史的行政区域データセットβ版』（CODH作成）
264 fields = None
265 records = []
266 for rec in _load_from_geoshape_city_csv():
267 # 「代表点」
268 refpt = sgeom.Point(rec.pop("longitude"), rec.pop("latitude"))
269 if not fields:
270 fields = list(sorted(rec.keys())) + ["refpoint"]
271 records.append([rec[k] for k in fields[:-1]] + [refpt])
272 #
273 cur = conn.cursor()
274 _do_insert(cur, "t_geoshape_city", fields, records)
275 #
276 # 駅データ.jp
277 # company
278 records = []
279 fields = [k for k, ty in coldefs_company]
280 for rec in _load_from_ekidata_csv("company20200619.csv", coldefs_company):
281 records.append([rec[k] for k in fields])
282 _do_insert(cur, "t_ekidata_company", fields, records)
283 # join
284 records = []
285 fields = [k for k, ty in coldefs_join]
286 for rec in _load_from_ekidata_csv("join20210312.csv", coldefs_join):
287 records.append([rec[k] for k in fields])
288 _do_insert(cur, "t_ekidata_join", fields, records)
289 # line
290 records = []
291 fields = None
292 for rec in _load_from_ekidata_csv("line20210312free.csv", coldefs_line):
293 refpt = sgeom.Point(rec.pop("lon"), rec.pop("lat"))
294 if not fields:
295 fields = list(sorted(rec.keys())) + ["refpoint"]
296 records.append([rec[k] for k in fields[:-1]] + [refpt])
297 _do_insert(cur, "t_ekidata_line", fields, records)
298 # station
299 records = []
300 fields = None
301 for rec in _load_from_ekidata_csv("station20210312free.csv", coldefs_station):
302 refpt = sgeom.Point(rec.pop("lon"), rec.pop("lat"))
303 if not fields:
304 fields = list(sorted(rec.keys())) + ["refpoint"]
305 records.append([rec[k] for k in fields[:-1]] + [refpt])
306 _do_insert(cur, "t_ekidata_station", fields, records)
307 #
308 conn.commit()
309
310
311 def _main(bb):
312 qrect = sgeom.box(*bb)
313 first = not os.path.exists(_DBNAME)
314 with spatialite.connect(_DBNAME, detect_types=sqlite3.PARSE_DECLTYPES) as conn:
315 if first:
316 _initialize_db(conn)
317 _fill_data(conn)
318 c = conn.cursor()
319 # 役立つものにするには、最低でも以下の検索をフレキシブルにせねばね。
320 c.execute("""
321 SELECT address, code, refpoint FROM t_geoshape_city
322 where ST_Within(refpoint, ?)""", (qrect,))
323 for found in c.fetchall():
324 print(found[:-1], found[-1].wkt)
325 c.execute("""
326 SELECT station_name, address, refpoint FROM t_ekidata_station
327 where ST_Within(refpoint, ?)""", (qrect,))
328 for found in c.fetchall():
329 print(found[:-1], found[-1].wkt)
330 c.execute("""
331 SELECT line_name, line_cd, refpoint FROM t_ekidata_line
332 where ST_Within(refpoint, ?)""", (qrect,))
333 for found in c.fetchall():
334 print(found[:-1], found[-1].wkt)
335
336
337 if __name__ == '__main__':
338 import json
339 import argparse
340 ap = argparse.ArgumentParser()
341 ap.add_argument("bb")
342 args = ap.parse_args()
343 _main(json.loads(args.bb))
雑にやってるし、まぁこういう風に混ぜ込んじゃって便利に感じるのは「オレさまだけ」なのであろうし、ほんとにこういうことをしだすと、openflight のデータだのの航空関係も入れたくなる、だの、キリがない。けどまぁ geoshape-city だけをやってると気付かない問題なんかにも出くわすわけで、似たことをしたい人には、多少は役に立つのではないかな、と思ってな。
ところで SQLiteStudio なのだけど、SELECT load_extension('mod_spatialite') がそのまま使えて、一応ちゃんと spatial を扱える SQLiteStudio となる:
ただ、テーブルビューで Geometry 型が解釈されることがないので、中身が空っぽにみえる。絵でやってるように AsWKT なりで明示的にお問合せないと中身がみえない。残念。
