2

I've converted a huge Shapefile into a SpatiaLite file because I need to make a query that uses a spatial index. My problem started on this question where I wanted to filter a specific area of a Shapefile once I give it a central point. To sum up, I already have the query that I need to use inside SpatiaLite to get the result that I want, and it's the following:

SELECT AsText(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));

So far everything is fine. Since I'm using an index inside Spatialite, I get the result of this query very fast (that doesn't happen working directly with the Shapefile). In this particular case, I'm filtering an area of my spatial data and printing only one column as the following:

POLYGON((339380.3182 7391671.753232, 339374.934121 7391681.671268, 339374.280746 7391682.874854, 339370.941839 7391689.025468, 339370.551171 7391689.745119, 339368.918082 7391692.75344, 339366.908313 7391691.723008, 339365.730444 7391691.1191, 339364.196167 7391690.332457, 339361.078168 7391695.89327, 339362.318072 7391696.583444, 339366.930247 7391699.069134, 339370.258422 7391700.862824, 339371.949654 7391701.774299, 339377.618093 7391691.15344, 339386.428107 7391675.083422, 339380.3182 7391671.753232))
POLYGON((339279.256409 7391651.056254, 339282.375504 7391645.321987, 339283.474543 7391643.301472, 339289.499822 7391632.224363, 339293.61682 7391624.655511, 339287.856062 7391621.52085, 339284.819854 7391619.849901, 339283.129676 7391622.964252, 339282.675141 7391623.801785, 339280.113922 7391628.521133, 339274.707503 7391638.483095, 339272.05454 7391643.371491, 339270.457411 7391646.314388, 339278.437978 7391650.623368, 339279.256409 7391651.056254))
POLYGON((339328.628779 7391643.652018, 339324.398045 7391641.363351, 339322.038033 7391645.683371, 339310.01731 7391667.687396, 339314.238657 7391670.032229, 339316.051101 7391666.709626, 339319.190366 7391660.954674, 339321.27978 7391657.124326, 339323.157066 7391653.682853, 339328.628779 7391643.652018))
POLYGON((339374.568446 7391654.673207, 339378.758524 7391652.424178, 339377.488583 7391650.053132, 339370.857976 7391637.673437, 339368.977998 7391634.163418, 339368.379575 7391633.046131, 339364.608089 7391635.123367, 339362.545075 7391636.197734, 339360.998084 7391637.003369, 339359.458082 7391638.01337, 339358.25808 7391639.053371, 339357.408079 7391639.873371, 339356.668077 7391643.633376, 339358.458079 7391645.873379, 339363.48115 7391648.621743, 339366.77357 7391650.420439, 339369.631314 7391651.981666, 339372.621277 7391653.611683, 339374.568446 7391654.673207))
POLYGON((339245.337594 7391672.994081, 339251.434035 7391660.501458, 339252.832571 7391657.635626, 339256.731889 7391649.645276, 339260.88073 7391641.143611, 339242.445386 7391631.18974, 339236.358771 7391643.483767, 339226.387344 7391663.624516, 339236.36559 7391668.558057, 339245.337594 7391672.994081))
POLYGON((339345.539182 7391631.400763, 339341.776631 7391622.903986, 339339.795756 7391618.473851, 339338.909817 7391616.492489, 339337.183689 7391612.63208, 339336.196467 7391610.424204, 339335.418578 7391608.684491, 339333.857078 7391609.521532, 339332.038896 7391610.496169, 339330.121376 7391626.79506, 339330.050935 7391627.393801, 339329.697556 7391630.397522, 339331.940096 7391631.651616, 339332.607729 7391631.967136, 339333.302411 7391632.220249, 339334.008241 7391632.441579, 339334.723724 7391632.630659, 339335.447345 7391632.787087, 339336.177572 7391632.910533, 339336.912858 7391633.000736, 339337.651648 7391633.057503, 339338.392376 7391633.080716, 339339.133474 7391633.070325, 339339.873373 7391633.026352, 339340.610506 7391632.94889, 339341.527311 7391632.801685, 339342.070241 7391632.694225, 339342.778782 7391632.483273, 339345.539182 7391631.400763))
POLYGON((339269.158185 7391645.612891, 339270.457411 7391646.314388, 339272.05454 7391643.371491, 339274.707503 7391638.483095, 339280.113922 7391628.521133, 339282.675141 7391623.801785, 339283.129676 7391622.964252, 339284.819854 7391619.849901, 339274.36854 7391614.091183, 339273.066244 7391619.062384, 339272.182833 7391622.434589, 339270.351498 7391629.425262, 339270.071456 7391630.494253, 339266.488471 7391644.171422, 339269.158185 7391645.612891))
POLYGON((339296.868015 7391660.383355, 339300.027383 7391654.713032, 339301.190816 7391652.624943, 339302.439346 7391650.384125, 339307.102077 7391642.015619, 339309.845443 7391637.091921, 339311.390204 7391634.319439, 339301.647443 7391629.024359, 339292.222049 7391645.835651, 339290.288803 7391649.283821, 339288.926974 7391651.712803, 339286.999041 7391655.151496, 339287.228276 7391655.272769, 339296.457124 7391660.155117, 339296.868015 7391660.383355))
POLYGON((339357.66806 7391728.533483, 339360.978064 7391730.173486, 339371.878888 7391710.113868, 339375.459772 7391703.52435, 339377.846385 7391699.132523, 339387.884077 7391680.661234, 339389.888111 7391676.973425, 339386.428107 7391675.083422, 339377.618093 7391691.15344, 339371.949654 7391701.774299, 339357.66806 7391728.533483))
POLYGON((339298.446288 7391719.026958, 339297.347986 7391721.293461, 339285.157265 7391745.884489, 339292.010873 7391749.300595, 339292.550974 7391749.663495, 339293.586875 7391750.167396, 339306.519997 7391723.117265, 339298.446288 7391719.026958))

Here I'm getting only the column with the polygons because that's the only thing I want inside this database.

My problem starts when I try to convert this data to GeoJSON, I know that I can use the ogr2ogr to convert CSV to GeoJSON, but that would require a specific syntax and only using ogr2ogr -f GeoJSON file.geojson file.csv wouldn't work. I also know that gdal has a VRT file that is supposed to be used with CSV when I want to convert them to GeoJSON, so I can specify the meaning of the columns as described on this question. However, right now I'm not even sure if the VRT file can do something about the SpatiaLite column syntax of representing polygons as POLYGON((339380.3182 7391671.753232, 339374.934121 7391681.671268, 339374.280746...)), or if I need to do it myself before using this VRT file. Also, I'm not sure if I'm overcomplicating the problem since the only thing I want is to get a GeoJSON file as a result from a query that I make on my SpatiLite database.

Is there any easier way of transforming the results of a SpatiaLite query into GeoJSON? Or transforming a SpatiaLite query on CSV, and then converting this CSV file to GeoJSON using a VRT file is the way to go?

Improve this question
2
  • ogr2ogr a.json file.sqlite
    – wingnut
    Commented May 7, 2021 at 5:31
  • 1
    For point data, if the X column starts with X, etc,ogr2ogr a.json input.csv -oo X_POSSIBLE_NAMES="X*" -oo Y_POSSIBLE_NAMES="Y*" -oo KEEP_GEOM_COLUMNS=NO should work. Shapefiles do in fact have indices, but you must create them. QGIS can do it, so can ogrinfo via ogrinfo -sql 'CREATE SPATIAL INDEX ON "layername"' shapefile.shp
    – wingnut
    Commented May 7, 2021 at 5:42

3 Answers 3

Reset to default
2

I think that converting data into SpatiaLite in this case may be a good idea, even shapefiles do support simple BBOX queries with spatial index if spatial index is created (the .qix or .sbn/.sbx files) and you can run your query with ogr2ogr also against the shapefile by using the SQLite SQL dialect. https://gdal.org/user/sql_sqlite_dialect.html.

You should get GeoJSON from ogr2ogr directly this way

ogr2ogr -f GeoJSON -sql "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))" output.json input.sqlite

You have something wrong in your query or in your data, or both. This is not correct

ST_Buffer(ST_GeomFromText('POINT ( 339298.446287788 7391719.02695803 )', 4326), 100)

EPSG:4326 coordinates are within range -90-> 90 and -180->180 degrees and a buffer that is 100 degrees wide means about 12000 kilometres. Check your data and what is the right coordinate system. You may get right results now but it is by chance.

Improve this answer
3
  • Thanks! The first command solved the issue for me... Actually, I know that I can create an index as a .qix file from a Shapefile. But apparently I can't use this index together with the SQLite SQL dialect. So in my particular case, the case was taking 6 seconds to finish. The original answer of the question I've mentioned talked about creating an index for the Shapefile. But the final conclusion was that it was not possible. That's why I was trying with SpatiaLite now.
    – raylight
    Commented May 7, 2021 at 15:09
  • About the query, I've tried creating a new Shapefile and SpatiaLite database changing the projection with -s_srs EPSG:31983 -t_srs EPSG:4326 and took out the 4326 from the ST_GeomFromText function, but for some reason it was not working... This is actually an issue I'm trying to solve. But using this query this way is the only way that is working so far (I can see the circle when I convert it to Shapefile and I open it with QGIS).
    – raylight
    Commented May 7, 2021 at 15:13
  • 1
    The SQLite dialect actually imports the data from the shapefile into a temporary SQLite database and the connection to the native .qix index is lost. It might be used with more simple use case like running ogr2ogr with -spat option, but maybe not even then.
    – user30184
    Commented May 7, 2021 at 19:26
1

Index the shapefile with:

ogrinfo -sql 'CREATE SPATIAL INDEX ON "layername"' shapefile.shp

To get the layer name, if you don't know it, use:

ogrinfo shapefilename.shp

eg,

1: 2019 Lines (Line String)

This means the first layer in the file is "2019 lines" and it's a linestring. The query is run on the layer, not the file. The file is considered the database.

Load point data as shapes from CSV with:

ogr2ogr a.json input.csv -oo X_POSSIBLE_NAMES="X*" -oo Y_POSSIBLE_NAMES="Y*" -oo KEEP_GEOM_COLUMNS=NO

assuming that the CSV file column names for east and north are something like X* and Y*

Improve this answer
6
  • Shapefile has always only one layer and the name of that is the same as the base filename. Thus the layer in shapefilename.shp is shapefilename but it is OK to check it with ogrinfo.
    – user30184
    Commented May 7, 2021 at 10:07
  • Thanks for pointing that out. Appending another shapefile renames the layer.
    – wingnut
    Commented May 7, 2021 at 11:09
  • Is your suggestion creating a .qix file and working directly with the Shapefile? I've tried that but I didn't manage to make it work, it seems that the index is not used when the SQLite SQL dialect is used... That was the first suggestion on this question but even though the query returns the result, it still takes too long to finish. While inside SpatiaLite the result is immediate.
    – raylight
    Commented May 7, 2021 at 15:17
  • Maybe don't use the Spatialite dialect?
    – wingnut
    Commented May 7, 2021 at 17:49
  • I'm not sure if I'm missing something... But if I don't use Spatialite dialect I receive an error saying Undefined function 'ST_GeomFromText' used. To sum up, if I don't use Spatialite dialect it won't recognize any of my spatial functions used in my query... When I'm working directly with a Spatialite database file instead of Shapefile that doesn't happen.
    – raylight
    Commented May 7, 2021 at 19:17
1

I just had the same problem, and for me using AsGeoJSON() instead of AsText() seemed the best option. See here: https://www.gaia-gis.it/gaia-sins/spatialite-sql-4.4.0.html

Improve this answer

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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