[GIS] Invalid join selectivity error in postgis 2.1.3

pointpolygonpostgiswithin

I have 2 tables, one for points and one for polygons. I want to find points falling within theese polygons, so I quered:

select p.id, l.id
from points p inner join polys l on st_contains(l.geom, p.geom);

and i got invalid join selectivity error. I know that was issue in previous version of postgis but as far as I know it should be fixed in 2.1.3.
In my tables I have about 400.000 points and 1.000 polygons – maybe too many?

Best Answer

I have had a similar problem. I believe this message is telling you that the projections of the two things that you want to join are incompatible.

In my case, it was a problem with the point data, which I had loaded into postgis using qgis. The translation of one to the other didn't work and so my points had nonsense geometry. This seems to be a more general problem with the way my version of qgis (2.4.0) reads in csv files and saves projections.

One solution is to try transforming your data in postgis so that both sides of the join have the same projection:

ALTER TABLE mytable
ALTER COLUMN geom TYPE geometry(Point,4283) 
USING ST_Transform(geom,4283);

Where 4283 is the projection.

Another solution is to re-create the points data in something else, for example in R:

### Libraries I always seem to use
library(sp);           ## Spatial object functions
library(rgdal);        ## Spatial data I/O utilities
library(rgeos);        ## GEOS spatial functions
library(raster);       ## Raster object functions
library(splancs);      ## Functions: n2dist (nearest neighbour)
library(reshape2)      ## Reshape datasets for simpler plotting
library(ggplot2);      ## ggplot used for majority of final plots
library(data.table)    ## used for big merges 


### read in your points:

mytable<- read.csv("filepath/filename.csv", header=TRUE, sep=",")


### coerce them into a spatial points data frame:

mytable<- SpatialPointsDataFramemytable[, c("x", "y")],
    data.frame(ID= mytable$id), ##an alternative if you don't have an id column: seq(1:nrow(mytable))
    proj4string=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")) 

### set the spatial points data frame to the projection you want:

mytable <- spTransform(mytable, CRS("+init=epsg:4283"))


### export them as a shape file:

writeOGR(mytable,
         "filepath/",
         "mytable",
         driver="ESRI Shapefile")
Related Question