I have need your help with an spatially SQL-query with spatialite. I would like to extract points from an sql database geospatially based on another polygon shape (only these points within the polygon box). I have found some approaches but the problem is that not only the expected 4 greenpoints are returned (see below) but instead the whole locations dataframe with all 7 points is returned. I would expect point IDs 1, 2, 3 and 4 to be returned - IDs 0, 5 and 6 are outside the bounding box.
It is plotted here (Fig: Locations within box), where the points are the locations and the box is the polygon. The expected result would be to only see all green points.
The locations inside the SQL data base is looking like this:
locations=
WKT ID
0 POINT (-1.2 11.0 0) 0
1 POINT (2.3 8.0 0) 1
2 POINT (1.4 8.5 0) 2
3 POINT (0.4 8.7 0) 3
4 POINT (3.7 9.6 0) 4
5 POINT (-0.4 7.2 0) 5
6 POINT (2.5 10.2 0) 6
The box inside the SQL data base is looking like this:
box=
index WKT
0 0 POLYGON ((0 6 0,5 6 0,5 10 0,0 10 0,0 6 0))
To select the locations within the box, I used the following query:
pd.read_sql_query("SELECT * FROM locations, box WHERE ST_Within(box.WKT, locations.WKT)", db)
which resulted in a join of both tables, but not in an selection of the locations table.
index WKT ID index WKT
0 0 POINT (-1.2 11.0 0) 0 0 POLYGON ((0 6 0,5 6 0,5 10 0,0 10 0,0 6 0)) (--> outside of box)
1 1 POINT (2.3 8.0 0) 1 0 POLYGON ((0 6 0,5 6 0,5 10 0,0 10 0,0 6 0))
2 2 POINT (1.4 8.5 0) 2 0 POLYGON ((0 6 0,5 6 0,5 10 0,0 10 0,0 6 0))
3 3 POINT (0.4 8.7 0) 3 0 POLYGON ((0 6 0,5 6 0,5 10 0,0 10 0,0 6 0))
4 4 POINT (3.7 9.6 0) 4 0 POLYGON ((0 6 0,5 6 0,5 10 0,0 10 0,0 6 0))
5 5 POINT (-0.4 7.2 0) 5 0 POLYGON ((0 6 0,5 6 0,5 10 0,0 10 0,0 6 0)) (--> outside of box)
6 6 POINT (2.5 10.2 0) 6 0 POLYGON ((0 6 0,5 6 0,5 10 0,0 10 0,0 6 0)) (--> outside of box)
Why am I getting a this result instead of an geospatially selection of the locations?
What I noticed is, that the geometry exportet as WKT(or WKB) does not have an assigend SRS after laoding from WKT (or WKB). Maybe this affects the bahavior?
.
.
.
.
The full script can be executed with:
import pandas as pd
from osgeo import ogr
import spatialite
# create a bounding box as WKT String and save as pd.DataFrame
box_df=pd.DataFrame(data=['POLYGON ((0 6 0,5 6 0,5 10 0,0 10 0,0 6 0))'], columns=["WKT"])
# generate some points, some within the bounding box and others outside
locations_df=pd.DataFrame(
data=[
'POINT (-1.2 11.0 0)',
'POINT (2.3 8.0 0)',
'POINT (1.4 8.5 0)',
'POINT (0.4 8.7 0)',
'POINT (3.7 9.6 0)',
'POINT (-0.4 7.2 0)',
'POINT (2.5 10.2 0)'
],
columns=['WKT'],
)
locations_df['ID']=range(len(locations_df))
# goal is to select only the yellow points within the box via an sql query
#sql script init
pathSQL = "testSQL2.sqlite"
db = spatialite.connect(pathSQL)
#write to db
locations_df.to_sql('locations', db)
box_df.to_sql('box', db)
#selection
pd.read_sql_query("SELECT * FROM locations, box WHERE ST_Within(locations.WKT, (SELECT WKT FROM box WHERE id=0))", db)
pd.read_sql_query("SELECT * FROM locations, box WHERE ST_Within(box.WKT, locations.WKT)", db)
#srs
print('SRS: ', ogr.CreateGeometryFromWkt(pd.read_sql_query("SELECT * from locations", db).WKT.iloc[0]).GetSpatialReference())