3

I have a SpatiaLite database using the projection EPSG:31983 and I can see it inside QGIS as the following:

enter image description here

When I use the following query I can successfully create a buffer of 100 meters on a point of this database:

SELECT GEOMETRY FROM lotes WHERE ST_Area(ST_Intersection(ST_Buffer(ST_GeomFromText('POINT ( 339298.446287788 7391719.02695803 )', 4326), 100), lotes.geometry)) > 0 AND lotes.rowid IN (SELECT rowid FROM SpatialIndex WHERE f_table_name = 'lotes' AND search_frame = ST_Buffer(ST_GeomFromText('POINT ( 339298.446287788 7391719.02695803 )', 4326), 100));

enter image description here

So far everything is fine...

My problem is that I'd like to use the projection EPSG:4326 instead of the EPSG:31983. So, I've converted my SpatiaLite database to the EPSG:4326 projection with the following command:

ogr2ogr -s_srs EPSG:31983 -t_srs EPSG:4326 -f SQLite -dsco SPATIALITE=YES lotes4326.sqlite lotes.sqlite

However, even though this conversion was successful, I don't manage to make my old query work properly after this projection conversion. The query either returns all my database or it returns nothing... I've tried the following:

SELECT GEOMETRY FROM lotes WHERE ST_Area(ST_Intersection(ST_Buffer(ST_GeomFromText('POINT ( -23.5770487595231 -46.57464981079102 )', 4326), 100), lotes.geometry)) > 0 AND lotes.rowid IN (SELECT rowid FROM SpatialIndex WHERE f_table_name = 'lotes' AND search_frame = ST_Buffer(ST_GeomFromText('POINT ( -23.5770487595231 -46.57464981079102 )', 4326), 100));

I've read this other issue that is very similar to mine where it's recommended to use points.geom::geography,...)::geometry or "just cast to geography". But I didn't manage to make it work on the context of my query. And since that question is about PostGIS and not SpatiaLite I'm not sure if the same solution applies to SpatiaLite... What am I missing? How can I use ST_Buffer in meters while using the projection EPSG:4326?

3
  • 1
    You have perfectly good data in metres why reproject to 4326 before buffering? – Ian Turton May 8 at 10:46
  • @IanTurton I want the coordinate system to be directly on 4326 just so my lat and long will be compatible with the lat and long that leaflet uses. In this example, before the reprojection the coordinates used was 339298.446287788 7391719.02695803. After the reprojection, the coordinates had a different number -23.5770487595231 -46.57464981079102 that's compatible with the coordinates leaflet captures. – raylight May 8 at 11:08
  • 1
    Leaflet (or any good WMS/WFS) is capable of reprojecting the data for you – Ian Turton May 8 at 15:50
4

First thing that makes fail all what you have tested is inversion of X and Y coordinates

'POINT( 339298.446287788 7391719.02695803) should become 'POINT ( -46.57464981079102 -23.5770487595231 )' and not 'POINT ( 23.5770487595231 -46.57464981079102)' like in your sample

Testable with

ogrinfo lotes.sqlite -dialect SQLite -sql "SELECT ST_Transform(ST_GeomFromText('POINT ( 339298.446287788 7391719.02695803 )', 31983), 4326)"

Or by going to https://epsg.io/transform#s_srs=31983&t_srs=4326&x=339298.4462878&y=7391719.0269580 (PS: longitude is X and latitude is Y)

Then to practically do the job, you have two approaches:

  • You can store your data in EPSG:31983 projection, use EPSG:4326 coordinates as an input, transform these input coordinates to local projection to create your buffer and return the geometry with EPSG:31983 projection transformed to EPSG:4326
ogrinfo lotes.sqlite -dialect SQLite -sql "SELECT ST_Transform(GEOMETRY, 4326) AS GEOMETRY FROM lotes WHERE ST_Area(ST_Intersection(ST_Buffer(ST_Transform(ST_GeomFromText('POINT ( -46.57464981079102 -23.5770487595231)', 4326), 31983), 100), lotes.geometry)) > 0 AND lotes.rowid IN (SELECT rowid FROM SpatialIndex WHERE f_table_name = 'lotes' AND search_frame = ST_Buffer(ST_Transform(ST_GeomFromText('POINT ( -46.57464981079102 -23.5770487595231 )', 4326), 31983), 100));"
  • You can store your data in EPSG:4326 "projection", transform these coordinates to local projection, do your buffer in the EPSG:31983 unit (meters), then transform back the buffer coordinates to EPSG:4326 to get "working" intersection and return the geometry with EPSG:4326
ogrinfo -so lotes4326.sqlite -dialect SQLite -sql "SELECT GEOMETRY FROM lotes WHERE ST_Area(ST_Intersection(ST_Transform(ST_Buffer(ST_Transform(ST_GeomFromText('POINT (-46.57464981079102 -23.5770487595231)', 4326), 31983), 100), 4326), lotes.geometry)) > 0 AND lotes.rowid IN (SELECT rowid FROM SpatialIndex WHERE f_table_name = 'lotes' AND search_frame = ST_Transform(ST_Buffer(ST_Transform(ST_GeomFromText('POINT (-46.57464981079102 -23.5770487595231)', 4326), 31983), 100), 4326));"

Edit:

As stated in @user30184 answer, these operations (transforming coordinates) to solve your issue are related to the fact buffer function is related to projection spatial reference unit from input geometry

2

The documentation of the SpatiaLite functions is in http://www.gaia-gis.it/gaia-sins/spatialite-sql-latest.html

Buffer

return a geometric object defined by buffering a distance around the geom, where dist is in the distance units for the Spatial Reference of geom.

It is not possible to directly define meters as unit when creating a buffer for EPSG:4326 geometry with SpatiaLite. You must transform the geometry into some metric CRS with ST_Transform first, create a buffer, and transform the result into EPSG:4326. You can chain all the functions into one SQL statement.

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.