导入Shapefile到数据库与我们程序中的其他版本差不多。
例如：
In [5]:
import sqlite3 as sqlite
db = sqlite.connect(':memory:')
db.enable_load_extension(True)
db.execute('SELECT load_extension("mod_spatialite.so")') # In Debian 9
cursor = db.cursor()
cursor.execute('SELECT InitSpatialMetaData();')
cursor.execute("DROP TABLE IF EXISTS gshhs")
cursor.execute("CREATE TABLE gshhs (" +
"id INTEGER PRIMARY KEY AUTOINCREMENT, " +
"level INTEGER)")
cursor.execute("CREATE INDEX gshhs_level on gshhs(level)")
cursor.execute("SELECT AddGeometryColumn('gshhs', 'geom', " +
"4326, 'POLYGON', 2)")
cursor.execute("SELECT CreateSpatialIndex('gshhs', 'geom')")
db.commit()
以下定义字符串模板以便后继的使用。注意里面的 WKT 。
In [6]:
sql_tpl = "INSERT INTO gshhs (level, geom) VALUES (2, GeomFromText('{0}', 4326))"
from osgeo import ogr
ogrfName = '/data/gdata/GSHHS_c.shp'
shapefile = ogr.Open(ogrfName)
layer = shapefile.GetLayer(0)
for i in range(layer.GetFeatureCount()):
feature = layer.GetFeature(i)
geometry = feature.GetGeometryRef()
wkt = geometry.ExportToWkt()
cursor.execute( sql_tpl.format(wkt))
db.commit()
In [7]:
import shapely.wkt
LONDON = 'POINT(-0.1263 51.4980)'
pt = shapely.wkt.loads(LONDON)
cursor.execute("SELECT id,level,AsText(geom) " +
"FROM gshhs WHERE id IN " +
"(SELECT pkid FROM idx_gshhs_geom" +
" WHERE xmin <= ? AND ? <= xmax" +
" AND ymin <= ? and ? <= ymax) " +
"AND Contains(geom, GeomFromText(?, 4326))",
(pt.x, pt.x, pt.y, pt.y, LONDON))
shoreline = None
for id,level,wkt in cursor:
#if level == 1:
shoreline = wkt
print(id, level, wkt)
15 2 POLYGON((-3.022472 58.643333, -4.438361 57.486194, -1.819167 57.610833, -2.53925 56.566639, -3.377528 56.380833, -2.595861 56.2675, -3.849139 56.117472, -2.137556 55.915806, -1.212556 54.580833, -0.077528 54.119944, 0 53.763879, 0.000028 53.76375, 0 53.763748, -0.780833 53.695833, -0.224194 53.650806, 0 53.543722, -0.007528 52.884972, 0 52.887507, 0.028278 52.897028, 1.301611 52.933722, 1.765806 52.478333, 0.932528 51.590806, 0.254167 51.465861, 1.426611 51.392889, 1.394194 51.153306, 0 50.787944, -0.005 50.787889, -3.476694 50.687056, -3.643361 50.218694, -5.715889 50.061667, -4.230917 51.186639, -3.056722 51.177056, -2.347583 51.7975, -3.403417 51.380389, -5.320833 51.860806, -3.9375 52.553333, -4.058361 52.922944, -4.769167 52.795, -4.199167 53.212472, -2.680861 53.354972, -3.107556 53.545861, -2.799278 54.239944, -3.151694 54.061194, -3.642528 54.508333, -3.052556 54.982472, -4.856722 54.86875, -4.856722 54.63125, -5.145833 54.854944, -4.884194 55.942472, -4.427528 55.903306, -4.75175 56.207083, -5.315889 55.851639, -5.805861 55.3025, -5.004167 56.712472, -6.229167 56.725833, -5.387528 57.108306, -5.814222 57.855, -5.072528 57.819139, -5.003361 58.627917, -3.022472 58.643333), (-4.334444 57.379444, -4.676389 57.138611, -4.335556 57.373056, -4.334444 57.379444))
In [8]:
f = open("xx_uk-shoreline.wkt", "w")
f.write(shoreline)
f.close()
In [9]:
cursor.execute("EXPLAIN QUERY PLAN " +
"SELECT id,level,AsText(geom) " +
"FROM gshhs WHERE id IN " +
"(SELECT pkid FROM idx_gshhs_geom" +
" WHERE xmin <= ? AND ? <= xmax" +
" AND ymin <= ? and ? <= ymax) " +
"AND Contains(geom, GeomFromText(?, 4326))",
(pt.x, pt.x, pt.y, pt.y, LONDON))
for row in cursor:
print (row)
(2, 0, 91, 'SEARCH gshhs USING INTEGER PRIMARY KEY (rowid=?)') (6, 0, 0, 'LIST SUBQUERY 1') (9, 6, 186, 'SCAN idx_gshhs_geom VIRTUAL TABLE INDEX 2:B0D1B2D3') (21, 6, 0, 'CREATE BLOOM FILTER')
运行上面的程序表明， SpatiaLite查询优化器将利用空间指数和数据表，并迅速识别边界盒的特征：
In [10]:
for row in cursor:
print(row)
(0, 0, 0, u'SEARCH TABLE gshhs USING INTEGER PRIMARY KEY (rowid=?)')
(0, 0, 0, u'EXECUTE LIST SUBQUERY 1')
(1, 0, 0, u'SCAN TABLE idx_gshhs_geom VIRTUAL TABLE INDEX 2:B0D1B2D3')
Out[10]:
(1, 0, 0, 'SCAN TABLE idx_gshhs_geom VIRTUAL TABLE INDEX 2:B0D1B2D3')