>>> from env_helper import info; info()
页面更新时间： 2024-07-23 22:05:13
运行环境：
Linux发行版本: Debian GNU/Linux 12 (bookworm)
操作系统内核: Linux-6.1.0-23-amd64-x86_64-with-glibc2.36
Python版本: 3.11.2
8.5. 导入Shapefile¶
导入Shapefile到数据库与我们程序中的其他版本差不多。
例如：
>>> import sqlite3 as sqlite
>>> db = sqlite.connect(':memory:')
>>> db.enable_load_extension(True)
>>> db.execute('SELECT load_extension("mod_spatialite.so.7")') # 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()
>>>
>>>
>>> sql_tpl = "INSERT INTO gshhs (level, geom) VALUES (2, GeomFromText('{0}', 4326))"
>>>
>>> from osgeo import ogr
>>> fName = '/gdata/GSHHS_c.shp'
>>> shapefile = ogr.Open(fName)
>>> 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()
8.5.1. 在表中进行空间查询查找¶
我们已经生成了数据库。 现在我们想从数据库中查找需要的多边形。下面是利用SpatiaLite来实现这一点的：
>>> import sqlite3 as sqlite
>>> db = sqlite.connect(':memory:')
>>> db.enable_load_extension(True)
>>> # db.execute('SELECT load_extension("libspatialite.so.5")') # In Debian 8
>>> db.execute('SELECT load_extension("mod_spatialite.so.7")') # 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()
>>>
>>>
>>> sql_tpl = "INSERT INTO gshhs (level, geom) VALUES (2, GeomFromText('{0}', 4326))"
>>>
>>> ogrfName = '/gdata/GSHHS_c.shp'
>>> shapefile = ogr.Open(fName)
>>> 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()
>>>
>>> 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
在缺省的情形下，
SpatiaLite进行查询的时候没有使用空间索引，因此我们在查找中不得不明确指出idx_gshhs_geom索引来优化查询过程。
然而要注意空间查询的机制，SpatiaLite不是利用Shapely来提取多边形来确保该点涵盖其中；
而是利用SpatiaLite的Contains()函数直接进行查找范围内整个多边形的检测。
8.5.2. 保存结果¶
通过前面的诸多步骤， 我们检索出了我们希望获取的结果。 我们将结果保存起来以便更进一步的使用。 保存的时候直接以纯文本的格式保存成 WKT 数据。
>>> f = open("uk-shoreline.wkt", "w")
>>> f.write(shoreline)
>>> f.close()
8.5.3. 优化¶
空间查询的具体实现是复杂的，然而理论上会产生一个快速而准确的答案。
已安装的SpatiaLite版本，可能无法使用SQLite命令行。因此，让我们利用Python的EXPLAIN QUERY PLAN命令：
>>> 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, 0, 'SEARCH gshhs USING INTEGER PRIMARY KEY (rowid=?)')
(6, 0, 0, 'LIST SUBQUERY 1')
(8, 6, 0, 'SCAN idx_gshhs_geom VIRTUAL TABLE INDEX 2:B0D1B2D3')
运行上面的程序表明， SpatiaLite查询优化器将利用空间指数和数据表，并迅速识别边界盒的特征：
>>> 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')
(1, 0, 0, 'SCAN TABLE idx_gshhs_geom VIRTUAL TABLE INDEX 2:B0D1B2D3')