Example efficient two pass Spatialite query to find if a point is inside a polygon

A spatial query can be expensive to run in terms of CPU processing power. In order to improve the query performance, typically, a first pass query is done to reduce the candidates first before running a second pass query to do actual spatial query. The following is an example using SpatiaLite to perform a two pass spatial query to find whether a point is contained within a polygon.

Point 1 is totally outside of any polygons; Point 2 is outside but within a bounding rectangle of a polygon; Point 3 is totally inside a polygon.

The two pass SpatiaLite query has the following syntax.
Note: the first approximate query uses the MBRCONTAINS function and the second finer query uses the ST_CONTAINS function.

SELECT ST_CONTAINS(polygonGeometry, geomfromtext('POINT(116.4688 -31.5502)', 4326))
FROM(
SELECT
polygonGeometry,
MBRCONTAINS(polygonGeometry, geomfromtext('POINT(116.4688 -31.5502)', 4326)) AS inMBR
FROM
(
SELECT geometry AS polygonGeometry FROM MyPolygon
)
WHERE inMBR = 1
)

Case 1: Totally outside 
For Point 1 (Totally outside), the query returns no rows since the point is outside of any bounding rectangles, as shown in the screenshot below.


Case 2: Outside of any polygons but within a bounding rectangle
For Point2, the query returns a result of 0 value, since the point is inside a bounding rectangle but outside of any polygon.

Case 3: Totally inside a polygon
For Point 3, the query returns a result of 1 since the point is inside a bounding rectangle and also inside a polygon.

Convert MultiLineString geometries to LineString using SpatiaLite SQL

Recently I made a mistake in creating a linear feature as a multi-line string instead of a plain line string geometry in SpatiaLite. I really did not want to digitize the features again so I decided to use SQL statements to convert the existing multi-line string geometries to line strings.

The following simple example shows how to do the conversion from a multi-line string table roads to a new line string table newroads. The basic steps involve (1) creating a new table, (2) populate the new table, (3) updating the new table with the line string geometry, (4) fixing the spatial metadata.

The table roads with the incorrect Multi-Line String geometry

Create a new table

1. In SpatiaLite, type in the following SQL statement to create a new empty table newroads from the roads table. Execute the statement.
   
   CREATE TABLE newroads AS SELECT * FROM roads LIMIT 0
   
   A blank table is created with the same structure as the roads table.
   
   

Populate the new table with the textual attribute values

1. In SpatiaLite, type in the following SQL statement to copy the rows (without the geometry field) from the original roads table into the new table newroads.
   
   INSERT INTO newroads(pkuid,name)SELECT pkuid,name FROM roads 
   
2. Execute the SQL statement.
   
   All the rows are copied into the new table newroads.
   
   A simple query of the new table newroads show that all fields are populated with the exception of the geometry field, which remains as null.
   
   
   
   

Update the geometry field

1. In SpatiaLite, type in the following SQL statement to form the SQL update statements to populate the geometry field.
   
   SELECT
   'UPDATE newroads SET Geometry='||'ST_GeomFromText('''||
   AsText(Geometry)||
   ''', 4326) WHERE pkuid='||
   pkuid||';' FROM roads
   
   Note: 4326 is the coordinate system reference used for the roads feature.
   
2. Execute the statement.
   
   The update SQL statements are generated.
   
   
3. In SpatiaLite GUI, click the grey area above and left of the resultant rows, as shown above.
   
   All the resultant rows are selected.
   
4. Mouse right click on the selected rows.
   
   A pop up menu appears.
   
   
5. Choose Copy.
   
   The resultant rows are copied to the OS Clipboard.
   
6. In a text editor, paste the rows into a new file, e.g. demorun.sql.
   
   
   
7. Change all occurrences of the string "MULTILINESTRING" to "LINESTRING".
   
8. Change all occurrences of the string "((" to "(".
   
9. Change all occurrences of the string "))" to ")".
   
   The statements should now look like the following example.
   
   UPDATE newroads SET Geometry=ST_GeomFromText('LINESTRING(-123.013367 49.328037, -122.980708 49.353917, -122.981324 49.386576, -122.975779 49.409992)', 4326) WHERE pkuid=1;
   
10. Save the changes.
    
11. In SpatiaLite, click the Execute SQL script button as shown below.
    
    
    
    The SQL Script Windows appears.
    
    
12. Select the newly created SQL script file, e.g. demorun.sql. Click Open.
    
    The geometry field in newroads is updated.
    
    
    
    

Fix the spatial metadata

1. In SpatiaLite GUI, mouse right click on the geometry field of the new table newroads.
   
   A pop up menu appears.
   
   
2. Choose Recover geometry column.
   
   The Recover geometry column dialog box appears.
   
3. In the SRID field, type in the coordinate system SRID value e.g. 4326.
   
4. In the Dims field, choose the appropriate dimension of the geometry field, e.g. XY.
   
5. In the Geometry Type list, choose LINESTRING.
   
6. Click OK.
   
   The geometry field is recovered.
   
   
   

Using SpatiaLite GUI to create a point geometry table from an ordinary table

Sometimes I have an ordinary table with numerical latitude and longitude columns and I need to construct a point geometry table from the records in the normal table, e.g. stations_raw. An example of such a table is shown in the screenshot below.

In the SpatiaLite GUI, enter the following SQL command to create a new table named stations:
CREATE TABLE stations AS 
SELECT
PK_UID,
code,
name,
ST_GeomFromText(
'POINT('||lng||' '||lat||')', 
4326
)
AS Geometry
FROM stations_raw

Note: 4326 in the example command is just the geographical coordinate system SRID of the data.

Execute the command.


The table stations is created.


The table is still a non-spatial table. In order to change it to a spatial table, the following steps need to be done.

1. In the SpatiaLite GUI, press mouse right click on the geometry column.
   
   
   
2. In the pop up menu, choose Recover geometry column,
   
   The Recover Geometry column dialog box appears.
   
   
3. In the SRID field, type in the data's SRID e.g. 4326 for geographical latitude and longitude data.
   
4. In the Dims field, choose the appropriate dimensions of the data, e.g. XY for 2-D.
   
5. In the Geometry field, choose the appropriate geometry of the data, e.g. POINT for point data.
   
6. Press OK.
   
   If the parameters are correct, the following message will appear.
   
   
   The GUI should show the table as a spatial table (with a globe).
   

Example SpatiaLite query to find neighbors of a polygon

A common spatial query is to find all the neighbor polygons touching a subject polygon. The screenshot below shows a SpatiaLite database of selected country feature polygons.

If you want to find all the neighboring countries of the country Canada, the following SpatiaLite SQL query can be used.
select b.name
from country a, country b
where touches(a.geometry, b.geometry)
and a.name = 'Canada'


The result 'USA' is returned.
 

SpatiaLite query to join point features by distance

Sometimes I have two related point features in a SpatiaLite database, one feature layer (nodes) contains the identifier I need, and the other feature layer (stations) contains the attributes I want. The only relationship I can make use of is the spatial proximity between the two feature layers. The following screenshots illustrate the problem.

nodes and stations feature layers are within a certain distance from each other

The nodes table contain the required id and geometry. 

The stations table contains the desired attributes

The problem can be solved by formulating a simple spatial join query in SpatiaLite. The following shows how to do the job (using QGIS to submit the query to the SpatiaLite database).

1. In QGIS, select Database | DB Manager | DB Manager.
   
   The DB Manager appears.
   
2. In the tree view, select the SpatiaLite database containing the two feature layers, e.g. vancity.sqlite.
3. Click SQL Window.
   
   The SQL  Window appears.
   
4. Type in the following SQL query.
   
   select b.id as id,
   a.code,
   a.name,
   y(transform(a.geometry, 4326)) as lat,
   x(transform(a.geometry,4326)) as lng
   from stations as a
   join nodes as b
   on (ST_DISTANCE(a.Geometry,b.Geometry) < 10)
   
   Note: it may be necessary to adjust the ST_DISTANCE function's maximum distance between both features so that it doesn't leave out or include features.
   
5. Click Execute.
   
   The join results appear.
   
6. Optional. If you want to display the results on the map, then the query should include the desired geometry field as shown below.
   
   select b.id as id, 
   a.code, 
   a.name, 
   y(transform(a.geometry, 4326)) as lat,
   x(transform(a.geometry,4326)) as lng,b.geometry
   from stations as a 
   join nodes as b 
   on (ST_DISTANCE(a.Geometry,b.Geometry) < 10)
   
   
   
   Toggle on Load as new layer. Choose the primary key e.g. id in the Column with unique integer values field. Choose the geometry field in the Geometry column field. Type in a Layer name (prefix) e.g. tmp. Click Load now!.
   
   The results layer is loaded in the map.
   

SpatiaLite SQL statement to find points in a line buffer

A typical GIS spatial query is to find all the point features in a buffer around a linear feature. This can be easily achieved using SpatiaLite database also. The following illustrates the workflow using two SpatiaLite feature tables - nodes and segments.

Using any SpatiaLite interface (graphical or command line interface), simply type in the following SQL statement to create a 50 meter buffer around the line segment feature using the spatial operator ST_Buffer:

SELECT ST_Buffer(Geometry, 50) AS buffer 
FROM segments
WHERE name='line1'


The screenshot below shows how the 50 meter buffer looks like on a map.

To identify all the point nodes within the buffer, we can wrap the buffer statement in a more complex SQL statement:

SELECT a.* FROM nodes AS a,
(
SELECT buffer FROM
(
SELECT ST_Buffer(Geometry, 50) AS buffer 
FROM segments
WHERE name='line1'
) AS b
)
WHERE ST_Contains(buffer, a.Geometry)

The screenshot shows how the results are like on a map.