0

I have this shapefile:

shapes.shp

in which I have imported in SpatiaLite using this command:

spatialite_tool --import --shapefile shapes --db-path db.sqlite --table shapes --charset UTF-8 --srid 27700 --coerce-2d -k

which provides me with a geometry column which seems valid:

SELECT Count(*), GeometryType("geometry"), Srid("geometry"), CoordDimension("geometry")
FROM "shapes"
GROUP BY 2, 3, 4

and the output:

Count(*)    GeometryType("geometry")  Srid("geometry")  CoordDimension("geometry")
----------  ------------------------  ----------------  --------------------------
509282      MULTIPOLYGON              27700             XY

I then transform the geometry to SRID 4326 using the following:

UPDATE shapes SET multi_geom = TRANSFORM(geometry, 4326);

I then run the follwing:

SELECT Count(*), GeometryType("multi_geom"), Srid("multi_geom"), CoordDimension("multi_geom") from shapes GROUP by  2, 3, 4;

Which proves me with this:

Count(*)    GeometryType("multi_geom")  Srid("multi_geom")  CoordDimension("multi_geom")
----------  --------------------------  ------------------  ----------------------------
509282      MULTIPOLYGON                4326                XY

All seems good.

Then I apply the Spatial Index via:

SELECT CreateSpatialIndex("shapes", "multi_geom");

So now I can query quite quickly any items that intersects a bounding box and it works well. The problem I now have is that I would like to display these golygons on a map using AsGeoJSON(), and though it works fine in most cases I do however have some very large polygons:

SELECT ST_NPoints(multi_geom) AS pointCount FROM shapes ORDER BY pointCount DESC limit 10;
pointCount
----------
460061    
303778    
248892    
241198    
238410    
238153    
205812    
193641    
192667    
178179 

Ultimately what I would like it to be able to use a form of "ST_Subdivide" function or something similar which will break up the polygons into different "rows" so that I am not sending 20MB of GeoJSON to a browser when a large polygon intersects with my bounding box. I can deal with some overlap, however some of the polygons span 10's of kilometres which in my eyes is wasted bandwidth and compute. I have been trying various different "methods" to try and achieve this to no avail, so I am certain I am missing something fundamental here.

An example of the sheer size of polygon could be seen running this query:

SELECT multi_geom
FROM shapes
WHERE ROWID IN (
        SELECT ROWID
        FROM SpatialIndex
        WHERE f_table_name = 'shapes'
            AND f_geometry_column = 'multi_geom'
            AND search_frame = BuildMbr(
                -2.022857666015625,
                52.109351808806345,
                -1.9955635070800781,
                52.11989325575632,
                4326
            )
            AND ST_Intersects(
        multi_geom,
                BuildMbr(
                    -2.022857666015625,
                    52.109351808806345,
                    -1.9955635070800781,
                    52.11989325575632,
                    4326
                )
            )
);

While in reality I would like the polygon broken up so I don't have to render it all.

1
  • We use a Question/Answer paradigm. Please do not place an Answer in the Question; instead create an Answer.
    – Vince
    Commented Feb 16, 2021 at 20:59

1 Answer 1

0

By the documentation SpatiaLite does have ST_Subdivide function.

http://www.gaia-gis.it/gaia-sins/spatialite-sql-latest.html

Divides geom into many parts until each part can be represented using no more than max_vertices.

If the optional argument max_vertices is not explicitly specified a limit of 128 vertices is implicitly assumed.

NULL will be returned on invalid arguments.

However, the function is not included in SpatiaLite 5.0.0-beta0 that I have. Perhaps it is included in the 5.0 release version.

EDIT

A straightforward option is to cut just the part or polygon that is covered by the bounding box of the map view. It can be done with function ST_Intersection.

return a geometric object that is the intersection of geometric objects geom1 and geom2

4
  • Subdividing and then evaluating against ST_Intersects causes issues due to the overlapping boundaries of polygons. I am also sure (though happy to be proven wrong) that subdividing against the large polygons will still return the entire polygon just chopped up and not saving myself an bandwidth or processing power. Unless I am missing something here I feel I need to split the broken up polygons on a per row basis to restrict wasted collections. Commented Feb 16, 2021 at 13:58
  • ST_Subdivide splits the whole polygon but you can create a new table from the splitted polygons as postgis.net/docs/ST_Subdivide.html and use that in the future. Or you could use ST_Intersection( polygon, bbox) on the server side and send just that the intersection as GeoJSON for saving bandwidth. Boundaries of subdivided polygons should not overlap but because of inaccuracy of floating point algebra it may be necessary to use ST_SnapToGrid. Can you give a concrete example of the problems with intersects?.
    – user30184
    Commented Feb 16, 2021 at 14:12
  • I am not sure if there is anything interesting for you in this document web.archive.org/web/20201026134241/http://www.gaia-gis.it/….
    – user30184
    Commented Feb 16, 2021 at 14:19
  • 1
    ST_Intersection was what I needed, it's taken my request from 17MB to 256KB. I'm pretty happy with this Commented Feb 16, 2021 at 15:45

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.