8

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.

3
  • Related stackoverflow.com/questions/9911650/…
    – Mapperz
    Mar 14, 2016 at 15:13
  • 1
    For the line, I'd start with WKT syntax, and then make a geometry from the text. There are examples of that in the Spatialite Cookbook. However you'll need to add the spatialite extension before spatialite functions will be available in sqlite. If that isn't enough help, leave a comment and I'll try to come up with something more descriptive. It might help to explain what the actual data is.
    – BradHards
    Mar 14, 2016 at 15:33
  • look at SpatiaLite and Python
    – gene
    Mar 14, 2016 at 15:50

1 Answer 1

18

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.

5
  • 2
    Just in case somebody gets this error like me in UBUNTU 18.04: Error: /usr/lib/x86_64-linux-gnu/libspatialite.so: undefined symbol: sqlite3_spatialite_init The library has to be loaded as: connR.execute('SELECT load_extension("mod_spatialite.so")') Jul 19, 2018 at 7:41
  • 1
    @RogerVeciana thanks for pointing this out. I added a brief clarification to the post above. Jul 20, 2018 at 9:47
  • 1
    For higher-level interfaces to SpatiaLite, you can turn to GeoPandas' PostGIS capabilities or SQLAlchemy with GeoAlchemy.
    – j08lue
    Dec 15, 2019 at 19:09
  • Not sure if this is a platform-specific thing, or if things have changed in the four years since this answer was posted, but for me on MacOS X with python 3.9 it no longer works - connection object has no attribute "enable_load_extension"
    – ibrewster
    Nov 11, 2021 at 17:53
  • It still works on Windows and Ubuntu with Python 3.9 in a conda environment. Nov 12, 2021 at 9:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.