How to use SpatiaLite functions in a Python script?

Question

I'm trying to find the intersection coordinates between a line (knowing its start and end points) and a polygon (SQlite file).

I've managed to make it work in QGIS by creating a shapefile for the line and then using DBmanager and the following SQL statement:

SELECT astext(st_intersection(polygon.geom, line.geom)) as XY
from polygon, line
WHERE st_intersects(polygon.geom, line.geom)

It works just fine but now I'm trying to make it work in a python script with the XY coordinates of the line and the polygon.sqlite file as inputs.

I'm trying to write the code and I've come to this so far:

#!/usr/bin/env python
# -*- coding: UTF-8 -*-

import gdal
import ogr
import sqlite3
import math
import shapefile
import pyspatialite



fichier=open("C:\\MyPath\\polygon.sqlite","r")
connR =sqlite3.connect("C:\\MyPath\\polygon.sqlite")

curR =connR.cursor()

#trying to add a table for the line:

curR.execute("""
               CREATE TABLE my_line(
                   id INTEGER PRIMARY KEY,
                   XMin INTEGER,
                   YMin INTEGER,
                   XMax INTEGER,
                   YMax INTEGER,
                   geom MULTILINESTRING
                )
""")
connR.commit

# user coordinates inputs

XA=input('X A = ');
YA=input('Y A = ');
XB=input('X B = ');
YB=input('Y B = ');

curR.execute("""INSERT INTO route VALUES (?,?,?,?,?,?)""", (1, XA, YA, XB, YB,"MULTILINESTRING"))
connR.commit()

#trying to use SpatiaLite functions:

curR.execute("""
        SELECT astext(st_intersection(polygon.geom, line.geom)) as XY
        from polygon, line
        WHERE st_intersects(polygon.geom, line.geom)
""")
connR.commit()

As you can see, my first problem is that I don't know how to create a "line" object in the table 'my_line'. My second problem is that I don't know how to use the SpatiaLite functions in a python script.

Answer 1

You could directly integrate your line into the database via a geometry column if you load spatialite as an extension of sqlite:

import sqlite3

connR = sqlite3.connect(':memory:')

connR.enable_load_extension(True)

#now we can load the extension
# depending on your OS and sqlite/spatialite version you might need to add 
# '.so' (Linux) or '.dll' (Windows) to the extension name

#mod_spatialite (recommended)
connR.execute('SELECT load_extension("mod_spatialite")')   
connR.execute('SELECT InitSpatialMetaData(1);')  

# libspatialite
connR.execute('SELECT load_extension("libspatialite")')
connR.execute('SELECT InitSpatialMetaData();')

curR = connR.cursor()


curR.execute('''CREATE TABLE my_line(id INTEGER PRIMARY KEY)''')
curR.execute('SELECT AddGeometryColumn("my_line","geom" , 4326, "LINESTRING", 2)')
connR.commit()

polygon_wkt = 'POLYGON ((11 50,11 51,12 51,12 50,11 50))'

XA = 11
YA = 52
XB = 12
YB = 49

line_wkt = 'LINESTRING({0} {1}, {2} {3})'.format(XA, YA, XB, YB)

curR.execute("""INSERT INTO my_line VALUES (?,GeomFromText(?, 4326))""", (1, line_wkt))

connR.commit()

cursor = connR.execute('''
    SELECT astext(st_intersection(geom, GeomFromText(?, 4326))) from my_line
    WHERE st_intersects(geom, GeomFromText(?, 4326))''', (polygon_wkt, polygon_wkt))

for item in cursor:
    print item

The code 4326 is the spatial reference ID and depicts a simple latitute-longitude reference system. See here. This of course would also work by storing the polygon in the database by inverting LINESTRING and POLYGON.