88

Different approach. Knowing that the pain is in ST_Intersection, and that true/false tests are fast, trying to minimize the amount of geometry passing through the intersection might speed things up. For example, parcels that are totally contained in a jurisdiction don't need to be clipped, but ST_Intersection will still probably go to the trouble of building ...


19

ST_Distance actually calculates the distance between all the pairs of points, so, as such, no index could be used. So your query will do a sequence scan and then choose those geometries that are less than the distance you specify away. You are looking for ST_DWithin, which does use an index. SELECT SUM(population) FROM points WHERE ST_DWithin(location, ...


13

For displaying purposes it is always good to use a spatial index. It will improve speed of both rendering and spatial queries. However, if you plan to update large quantities of objects, it might be wise to remove the spatial index during the update. Otherwise the update process will become significantly slower, because with every update the spatial index ...


12

This is an exciting question! How big is the raster you want to query? WKTRaster is stored in the database as a BLOB. In order to find the value at a specific point, from a known (x_0, y_0) corner coordinate row/column indices (i, j) are computed using (dx, dy) steps and rotation. With (i, j) known, the ST_Value() function can access the actual data at the ...


12

That's the gist of it. The R-tree allows you to make a very fast first pass and gives you a set of results that will have "false positives" (bounding boxes may intersect when the geometries precisely do not). Then you go over the set of candidates (fetching them from the shapefile by their index) and do a mathematically precise intersection test using, e.g., ...


12

# assume a list of feature ids returned from index and a QgsVectorLayer 'lyr' fids = [1, 2, 4] request = QgsFeatureRequest() request.setFilterFids(fids) features = lyr.getFeatures(request) # can now iterate and do fun stuff: for feature in features: print feature.id(), feature 1 <qgis._core.QgsFeature object at ...


11

If you want to create a spatial index for all features in your layer, you could try using the following: layer.dataProvider().createSpatialIndex() This will create the .qix spatial index file. Edit: As mentioned by Matthias Kuhn in his comments below, the QgsSpatialIndex class is only used temporarily as seen in the blog you linked to; whereas the ...


9

You've almost got it, but you've made a small error. You need to use the intersection method on the spatial index, rather than passing the index to the intersection method on the buffered point. Once you've found a list of features where the bounding boxes overlap, then you need to check if your buffered point actually intersects the geometries. import ...


9

It's possible in QGIS3. You can click "Create Spatial index" from the Layer Properties Dialog:


8

If you want to batch create indexes on geometry columns, you could try this plpgsql function I've just knocked up: CREATE OR REPLACE FUNCTION BatchIndex(sn text, tn text, cn text) RETURNS void AS $$ DECLARE i_exists integer; DECLARE idxname text; BEGIN idxname := 'idx_' || tn || '_' || cn; select into i_exists count(*) from pg_class where relname = ...


8

Here are two non-Spatialite solutions. For Rtree with Python, try this example with 13000 points -- should take a few seconds: from random import randrange from rtree import index from math import sqrt # Create a 3D index p = index.Property() p.dimension = 3 idx3d = index.Index(properties=p) # Make and index random data coords = [] for id in range(13000):...


8

There's two things going on here: the GIST API in PostgreSQL and the bindings of types to that API for the purposes of building an R-Tree. PostGIS necessarily uses the PostgreSQL GIST API. That's what it's for. That way we don't have to worry about transaction management or writing things to disk or all the other messy important things involved in ...


8

The description from Boundless is very simplified. It somehow misses the first important step: The index will be able to filter most points/bounding boxes before running any calculation because of their location in the R-Tree. R-Trees break up data into rectangles, and sub-rectangles, and sub-sub rectangles, etc. The index will look in which rectangle the ...


7

You're correct that your problem is that in order to sort, then limit to the closest, you have to generate the distance to each of the 100k points, which is extremely time consuming. Luckily, there's help. First, you want to make sure that your geometry column is indexed, if it's not, then you need to index it, or this won't work. Then you want to use the ...


7

In a blog post on the subject, Nathan Woodrow provides the following code: layer = qgis.utils.iface.activeLayer() # Select all features along with their attributes allAttrs = layer.pendingAllAttributesList() layer.select(allAttrs) # Get all the features to start allfeatures = {feature.id(): feature for (feature) in layer} def noindex(): for feature in ...


6

A standard multicolumn b-tree index on the two columns is probably the most effective solution provided that: both columns are used together in the where expression. A multicolumn index will not be used when only the second column is present A multicolumn index is more efficient than two indexes on both columns because it will use less storage and will be ...


6

An index is used to filter a query using the "WHERE" part of the statement. A GiST index is not used by ST_Union. Without a "WHERE" part, then no filtering (or index) is required to return the result, and the query just chugs through all the rows in the query. As described in the manual, not all functions make use of indicies, for example ST_Distance and ...


6

Yes. It does as it appears from looking at the source code for the Spatialite Data Provider. The QgsSpatiaLiteFeatureIterator class is the one that supplies the features to the map upon sending a rectangle extent. You can just search for 'spatialIndex' in that class to see they actually use the index if available.


6

I have found that rearranging the query so that the sub-query is at the same level as the initial select, essentially a Cartesian product, but then using the where clause to restrict the records read, will cause the indexes to be used and avoid a full table scan. SELECT * FROM osm_addr2 AS addr, (SELECT geometry FROM osm_addr2 WHERE osm_id=-332537) ...


6

PostgreSQL doesn't use indexes for functions, it uses indexes for operators only. What happens is function inlining. ST_INTERSECTS is defined as: CREATE OR REPLACE FUNCTION ST_Intersects(geom1 geometry, geom2 geometry) RETURNS boolean AS 'SELECT $1 && $2 AND _ST_Intersects($1,$2)' LANGUAGE 'sql' IMMUTABLE; And so the query gets rewritten to use the ...


5

Two aspects that I found sped up my PostGIS raster calculations, were using integer values in the raster, and using multi-band rasters where possible. In this case, can the DN value be stored as integers, if this is not already being done? The other thought (and I'm not certain it is relevant here) is to use multi-band rasters. For example, if you are ...


5

Something must be wrong with your mysql installation or the .ini settings. Just tested a geospatial index on my old mac (10.6.8 / MySQL 5.2). That configuration is similar to yours and I tested the big geodata dump (9 million records). I did this query: SET @radius = 30; SET @center = GeomFromText('POINT(51.51359 7.465425)'); SET @r = @radius/69.1; SET @...


5

On large databases or a database with may changes it can be very important to have spatial indexes in place and updated regularly. (Keeping it simple here) For example for Oracle Spatial indexing capabilities into the Oracle database engine is a key feature of the Spatial product. A spatial index, like any other index, provides a mechanism to limit ...


5

For comparisons, look at More Efficient Spatial join in Python without QGIS, ArcGIS, PostGIS, etc. The solution presented use the Python modules Fiona, Shapely and rtree (Spatial Index). With PyQGIS and the same example two layers, point and polygon: 1) Without a spatial index: polygons = [feature for feature in polygon.getFeatures()] points = [feature ...


5

FME should support QIX index files for Shape in FME 2019.1 For now, if you want FME to create a QIX file for QGIS then use a FeatureWriter followed by ogrinfo in a SystemCaller.


4

What would happen if you omit the "st_multi(st_intersection(a.geom,b.geom))" part? Doesn't the below query mean the same thing without it? I ran it on the data you provided. INSERT INTO parcel_jurisdictions(parcel_gid,jurisdiction_gid,isect_geom) SELECT a.orig_gid parcel_gid, b.orig_gid jurisdiction_gid, a.geom FROM valid_parcels a, ...


4

Issue As suggested by @Kelso, you can check the Execution Plan to see where your issue is. Your query is relatively complex, so it's a good idea. If we mouse over the Filter Cost which is showing 97%, we can see that it's 'stuck' on the STDistance Query. This means we need to find another way around it. Suggested Solution In short, my suggestion is to ...


4

I would not use a conventional point->polygon process because that expects your points to define the boundary of a polygon and it doesn't sound like yours do. It sounds like yours are hotspots that are somehow related. However, there are lots of ways to create polygons for this sort of situation depending on what is sensitive in your data. Here's a few ...


4

Effectively forcing the planner to do the thing you want might help. In this case, sub-setting the polygon table prior to executing the spatial join with the points table. You might be able to outwit the planner using "WITH" syntax: WITH polys AS ( SELECT * FROM area WHERE area.id in IN(28,29,30,31,32,33,25,26,27,18,19,20,21,12,13,14,15,16,17,34,35,1,2,...


4

Your questions are answered in Rtree's docs, which come with the source and are online at http://toblerity.github.com/rtree/tutorial.html. You can make your shapefile indexing code much simpler: for i, shape in enumerate(shapes): idx.insert(i, shape.bbox) Use Python's enumerate() function for this kind of thing whenever possible. It's tidy and ...


Only top voted, non community-wiki answers of a minimum length are eligible