[PYTHON] Try using SpatiaLite to store spatial information such as maps in SQLite

Overview

SpatiaLite is an extension of SQLite that can store spatial information such as maps. Geometry information written in Shp files and GeoJson can be stored and used in the database.

http://www.gaia-gis.it/gaia-sins/index.html

Distribution file description

SpatiaLite You can operate the database from the command line.

** Binary distribution ** For Windows, you can get the binaries of any platform from the following. http://www.gaia-gis.it/gaia-sins/index.html

mod_spatialite.dll/mod_spatialite.so SQLite extension module. SpatiaLite functions can be used by loading a DLL or so file using load_extension on SQLite.

** Binary distribution ** For Windows, you can get the binaries of any platform from the following. http://www.gaia-gis.it/gaia-sins/index.html

** Source code distribution ** You can create it by downloading and building the libspatialite source code from the following. https://www.gaia-gis.it/fossil/libspatialite/index

spatialite_gui You can use the functions of spatialite in the GUI.

For Geometry type columns such as POLYGON, you can check the contents as an image.

** Binary distribution ** For Windows, you can get the binaries of any platform from the following. http://www.gaia-gis.it/gaia-sins/index.html

build libspatialite

If you can't get mod_spatialite in binary, you need to build libspatialite.

The following libraries are required to build libspatialite-4.2.0. proj-4.8.0 proj.4 is a mapping projection library. http://trac.osgeo.org/proj/

geos-3.4.2 Geometry Engin is a port of Java Topology Suite to C ++, an API of spatial predicates and features for processing Geometry. http://trac.osgeo.org/geos/

After installing the above, execute the following command in the folder containing the libsatialite source.

./configure --disable-freexl
make
make install

--disable-free xl disables the function of linking with Excel. If you enable this, you will need freexl. https://www.gaia-gis.it/fossil/freexl/index

If you get an error on FreeBSD

You may get the following error on FreeBSD

/usr/bin/ld: cannot find -ldl

In FreeBSD, dlopen is a standard library, so -ldl is not needed. Search for -ldl in src / Makefile, delete it, and then try make again.

reference: http://sourceforge.net/p/idjc/discussion/458834/thread/246f0841/

Use mod_spatialite from SQLite.

As mentioned earlier, you can use spatialite with load_extension, but there are some caveats.

-SQLite version must be 3.7.17 or later.

· All dependent DLLs must be located in a location visible to SQLite. Specify the folder containing the related DLL with the environment variable PATH.

-SQLite and mod_spatialite processes must be built with the same bits. If SQLite is running on 32bit, mod_spatialite needs to be built on 32bit, and if SQLite is running on 64bit, mod_spatialite needs to be built on 64bit. If you need 64-bit binary SQLite on Windows, you can create it by following the article below. http://qiita.com/akaneko3/items/0e99c3c1366dfbad006f

-By loading the DLL or so file using load_extension, the spatialite function can be used thereafter.

SELECT load_extension('/usr/local/lib/mod_spatialite.so');

-When creating a database from SQLite3, the metadata managed by SpatiaLite has not been created. Execute the following SQL to create this metadata.

SELECT InitSpatialMetaData();

tutorial

Here, we will actually store the administrative area of the national land numerical information in spatialite and use it.

Acquire administrative areas nationwide with national land numerical information.

First, download the administrative districts nationwide from the page below. http://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-N03.html

Create an empty SQLite file.

Start spatialite_gui and create an empty SQLite database file.

Press the Creating a New (empty) SQLite DB icon and enter a file name.

If created successfully, the path to the created database will be displayed.

Up to this point, it is no different from a normal database.

Import shp file

Import the shp file of national land numerical information into the database.

Press the [Load Shape file] icon and select "N03-14_140401.shp" in the national land numerical information.

At this time, you will be asked for the encoding, so select "CP932".

Upon successful completion, a table called "N03-14_140401" will be created.

Check the table created from the shp file

Here, we will check the table created from the shp file. To see the columns in the created table, select the table, right-click and select Show columns.

N03_001 to N03_004 are prefectures, wards, cities, and N03_007 is a TEXT type column in which the code is stored. Geometry has the shape of the administrative area stored as a Blob and is the type specified by spatiaLite called POLYGON.

You can see the actual SQL statement used to create this table by right-clicking on the table and selecting Show CREATE statement.

CREATE TABLE "N03-14_140401" (
"PK_UID" INTEGER PRIMARY KEY AUTOINCREMENT,
"N03_001" TEXT,
"N03_002" TEXT,
"N03_003" TEXT,
"N03_004" TEXT,
"N03_007" TEXT, 
"Geometry" POLYGON)

Let's actually look at the contents of the table. To do this, execute the SQL statement.

SELECT * FROM "N03-14_140401" LIMIT 10;

When executed, the following information can be obtained.

Data is expressed in text from PK_UID to N03_007. However, since Geometry is BLOB data, people cannot see the contents.

There are two ways to check this content. One is to use a BLOB explorer. Select any cell in the Geometry column, right click and select "BLOB explore".

BLOB explore can represent data in various formats such as binary, image, SVG and GeoJSON.

One way to avoid using BLOB explore is to use SQL to display Geometry columns in any format.

The example below displays the Geometry column in text format.

--Since the data is long, display up to the first 100 characters
SELECT substr(ASTEXT(Geometry), 1, 100) FROM "N03-14_140401" LIMIT 10;

I was able to display POLYGON information as characters in ASTEXT. It can be output in any format. For example, if it is ASGEOJSON, it will be output as GeoJSON.

SELECT substr(ASGEOJSON(Geometry), 1, 100) FROM "N03-14_140401" LIMIT 10;

The functions provided by SpatiaLite of ASTEXT and ASGEOJSON are described in the "Spatial SQL functions reference guide" below.

http://www.gaia-gis.it/gaia-sins/spatialite-sql-4.2.0.html

Try to create a table with coordinate information by yourself

In the previous example, I created a table from a shp file, but this time I will create a table with coordinate information by myself.

In the table below, create a table with coordinates of feature information such as Tokyo Tower.

First, create only columns that do not contain feature information.

CREATE TABLE "places" (
"PK_UID" INTEGER PRIMARY KEY AUTOINCREMENT,
"name" TEXT)

After executing SQL, refresh will add the table.

Next, create a column representing feature information using AddGeometryColumn ().

Select AddGeometryColumn ('places', 'Geometry', 0, 'POINT', 'XY')

After executing SQL, refreshing will add a column and add a globe icon to the table.

Once you've created the table, it's time to enter the data. In the following example, the coordinates of the Tokyo Metropolitan Government are stored in the table.

INSERT INTO places (name, Geometry) VALUES(
  'Tokyo Metropolitan Government',
  GeomFromText('POINT(139.692101 35.689634 )')
)

GeomFromText is a function that changes a string to Geometry. This time I converted from text format, but you can create it from GeoJSON or Kml using GeoJSON, GeoMFromKml, etc.

If you create it all at once as shown below, it will work like that at first glance, but it will not work properly when creating an RTree index because the metadata is not associated with it.

CREATE TABLE "places" (
"PK_UID" INTEGER PRIMARY KEY AUTOINCREMENT,
"name" TEXT,
"Geometry" POINT)

Search for coordinate information

Here, we will search for coordinate information. Try to find out which administrative division contains the coordinates specified in the places table.

SELECT
  name, "N03-14_140401".*
FROM
 "N03-14_140401"
INNER JOIN places ON
 Contains("N03-14_140401".Geometry, places.Geometry)

The Contains function checks if the range of the first argument includes the area of the second argument. The result is as follows:

Since the Tokyo Metropolitan Government Building is located in Shinjuku Ward, the results can be said to be as expected. In the previous example, the data stored in the table was searched, but it is also possible to specify it directly in the WHERE area.

SELECT
  *
FROM
 "N03-14_140401"
WHERE
  Contains("N03-14_140401".Geometry,GeomFromText('POINT(139.692101 35.689634 )'))

MBR and rough search

The previous search was an accurate search for each coordinate, but you may want to get a rough position quickly.

In this case, use Minimum Bounding Rectangle --MBR. MBR represents the approximate size of GEOMETRY.

Consider the following complex shapes as red border shapes.

This allows for an inaccurate but fast search. Use Envelope () to get the MBR from GEOMETRY.

SELECT
   ASTEXT(Envelope(Geometry))
FROM
 "N03-14_140401"
LIMIT 10

The result is a simple shape like this:

Now, let's search for the coordinates including the same Tokyo Metropolitan Government Building using the MBR.

SELECT
  *
FROM
 "N03-14_140401"
WHERE
  MBRContains("N03-14_140401".Geometry,GeomFromText('POINT(139.692101 35.689634 )'))

Instead of improving the search speed, unnecessary data has also been extracted as shown below.

Search using RTree index

Normally, when dealing with DB, the search speed is improved by using the index. Unfortunately, you can't give Geometry an index that you can give to a letter or number. However, SpatiaLite leverages SQLite's RTree to support RTree indexes. Please refer to the following for the algorithm of RTree index.

Wonderful R*Tree Spatial Index http://www.gaia-gis.it/gaia-sins/spatialite-cookbook/html/rtree.html

As an overview of RTree, Geometry is made into a rectangle using MBR, and a tree is created according to the intersection of the rectangles to make it easier to search. It is called RTree index because it creates a tree using rectangle.

Creating and using RTree indexes

Now let's actually use the RTree index. Use the CreateSpatialIndex function to create an RTree index. Specify the table name in the first argument and the column name in the second argument.

SELECT CreateSpatialIndex("N03-14_140401", "Geometry") ;

If it was created successfully, 1 will be returned and you can see that the SpatialIndex shows the new table after refreshing.

By creating an RTree index, the following four tables will be created.

・ Idx_N03-14_140401_Geometry ・ Idx_N03-14_140401_Geometry_node ・ Idx_N03-14_140401_Geometry_parent ・ Idx_N03-14_140401_Geometry_rowid

Except for idx_N03-14_140401_Geometry, it is a table used internally and we do not expect the user to operate it directly.

idx_N03-14_140401_Geometry is a VIRTUAL TABLE, which actually uses an internal table to return the result.

Now, let's check the contents of idx_N03-14_140401_Geometry.

SELECT * FROM "idx_N03-14_140401_Geometry" LIMIT 10;

The MRB of each record is stored in this way.

By using this table, users can filter data from a large amount of data and obtain the result.

SELECT 
  * 
FROM "N03-14_140401"
WHERE ROWID IN(
  SELECT
    pkid
  FROM
   "idx_N03-14_140401_Geometry"
  WHERE
    MBRContains(
      BuildMBR(xmin,ymin,xmax,ymax),
      GeomFromText('POINT(139.692101 35.689634 )')
  )
) AND Contains(Geometry, 
  GeomFromText('POINT(139.692101 35.689634 )')
)

I think I was able to get accurate data faster than a normal search. The larger the number of data, the greater the difference.

In this way, unlike a normal index, the RTree index is not applied just by executing SQL on the real table, but it has the role of filtering data by subquery or JOIN using the created VIRTUAL TABLE.

By managing the index in a separate table, you may be wondering if there may be a discrepancy when the real table is updated. However, since the trigger was created when CreateSpatialIndex was executed, the user can synchronize the index and the table without any modification.

To see this trigger, run the following SQL:

select name  from sqlite_master where type = 'trigger' and tbl_name='N03-14_140401';

This time, you can see that the following trigger is automatically created.

・ Ggi_N03-14_140401_Geometry ・ Ggu_N03-14_140401_Geometry ・ Gii_N03-14_140401_Geometry ・ Giu_N03-14_140401_Geometry ・ Gid_N03-14_140401_Geometry

Delete RTree index

Use DisableSpatialIndex to delete the RTree index you created.

SELECT DisableSpatialIndex("N03-14_140401", "Geometry") ;

You can see that the Spatial Index has been deleted by executing refresh.

Use from python

When using SpatiaLite from python, it's easy to use if SQLite is properly available.

You can find out which version of Sqlite3 Python is using by following these steps:

import sqlite3
print (sqlite3.sqlite_version_info)

If it does not meet the version, please update SQLITE. For Windows, you should rewrite the DLL in the following folder.

C:\Python27\DLLs\sqlite3.dll

Below is the code in Python.

#!/usr/bin/python
# -*- coding: utf-8 -*-
import sqlite3
import os

# mod_Add the folder with spatialite to your PATH
os.environ["PATH"] = os.environ["PATH"] + ';C:\\tool\\spatialite\\mod_spatialite-4.2.0-win-x86'
cnn = sqlite3.connect('database/gyouseikuiki.sqlite')

# mod_Loading spatialite
cnn.enable_load_extension(True)
cnn.execute("SELECT load_extension('./mod_spatialite-4.2.0-win-x86/mod_spatialite.dll');")
sql = """
SELECT
  N03_001,
  N03_002,
  N03_003,
  N03_004,
  N03_007, 
  AsGeoJson(Geometry)
FROM
 "N03-14_140401"
WHERE
  MBRContains("N03-14_140401".Geometry,GeomFromText('POINT(139.692101 35.689634 )'))
"""

ret = cnn.execute(sql)
for r in ret:
  print('----------------------------')
  print(r[0].encode('utf_8') )
  print(r[1].encode('utf_8') )
  print(r[2].encode('utf_8') )
  print(r[3].encode('utf_8') )
  print(r[4].encode('utf_8') )
  print(r[5].encode('utf_8') )

Pass the environment variable through the path where mod_spatialite is, allow_load_extension to use the extension DLL, and then load mod_spatialite.

After that, you can execute SQL as usual.

If an error occurs with enable_load_extension

The following error may occur when using enable_load_extension.

"'sqlite3.Connection' object has no attribute 'enable_load_extension'"

The reason is that Python, which is pre-installed on some operating systems such as MaxOS and Debian, disables the extension library loading feature at compile time.

https://docs.python.org/2/library/sqlite3.html#f1

The only way to deal with this is to reinstall Python. During the installation, you will need to manually modify setup.py in the Python source code folder.

Search for SQLITE_OMIT_LOAD_EXTENSION, comment out that line, and run ./configure, make, make install.

If you get an error with make install, try "make -i alt install" if you get an error during make install.

http://python.g.hatena.ne.jp/nelnal_programing/20101026/1288084718

How to use Peewee

Please refer to the following for how to use spatialite with the Peewee library which is an ORM.

http://qiita.com/mima_ita/items/9d4e1d0afac1865acdbb#spatialitesql%E3%81%B8%E3%81%AE%E6%8E%A5%E7%B6%9A%E6%96%B9%E6%B3%95

Use from php

In the case of Windows + PHP, it is not easy to use. See below for details.

I have a hard time using SpatiaLite when trying to store spatial information such as maps in the DB with PHP http://qiita.com/mima_ita/items/1a90de89f0a194ba843e

Other Tips

Perform geodetic system conversion using SpatiaLite

You can use the Transform function to transform the datum. In the example below, JGD2000 / plane orthogonal coordinate system 6 is converted to a world geodetic system.

select AsText(Transform(GeomFromText('POINT(-4408.916645 -108767.765479)', 2448), 4326))

reference

SpatiaLite 4.2.0 SQL functions reference list http://www.gaia-gis.it/gaia-sins/spatialite-sql-4.2.0.html

The SpatiaLite Cookbook http://www.gaia-gis.it/gaia-sins/spatialite-cookbook/index.html

** A quick tutorial to SpatiaLite --a Spatial extension for SQLite (old but useful) ** http://www.gaia-gis.it/gaia-sins/spatialite-tutorial-2.3.1.html