Why am I getting SRID error when doing spatial join in PostGIS?

47 views Asked by At

I have two tables that I wish to join:

  • sd_lat_long: This table contains UK locations and their respective latitudes and longitudes
  • sd_geo_uk_counties: This table contains UK counties and their respective geometries

I wish to create a view which finds which UK locations intersect with which counties.

Here is a script that I have come up with:

--Define list of uk locations and their lat, longs
WITH geo AS ( 
    SELECT 
        uk_locations,
        ST_SetSRID(ST_Point(latitude ,longitude),4326) AS lat_long --set to SRID 4326
    FROM
        sd_lat_long 
),
--Defining the shapefile and transforming the geometry
shp AS (
SELECT 
    county,
    ST_Transform(geometry, 4326) AS geometry --Transformed to SRID 4326 from 27700
FROM 
    sd_geo_uk_counties
WHERE 1=1
)
--Join on intersection of lat_long and the county geometry to assign a county to the UK locations
SELECT
    geo.uk_locations,
    shp.county
FROM 
    geo
LEFT JOIN
     shp
ON
    ST_Intersects(shp.geometry, geo.lat_long) 

For the shp CTE, I used this post to help: Why are UK county geometries showing up in the Gulf of Guinea when using PostGIS?

The problem that I'm facing is that I keep getting an error that appears to be related to the SRID of the geometries that I am using in my script. See error below. enter image description here

I firstly had a look at the SRIDs for the geometry fields in the geo and shp CTEs for the script above.

GEO: The lat_long field looked to be fine - i.e. SRID = 4326 Test script:

WITH geo AS ( 
    SELECT 
        uk_locations,
        ST_SetSRID(ST_Point(latitude ,longitude),4326) AS lat_long --set to SRID 4326
    FROM
        sd_lat_long 
),  
SELECT
    uk_locations AS geo_field,
    ST_SRID(lat_long) AS srid,
    'geo' AS cte
FROM
    geo
GROUP BY 
1,2, 3
HAVING
srid <> 4326

This returned no errors or entries with SRID <> 4326

enter image description here

SHP: Test script:

WITH shp AS (
SELECT 
    ctyua23nm AS county,
    ST_Transform(geometry, 4326) AS geometry
FROM 
    sd_geo_uk_counties
WHERE 1=1
)
SELECT
    county AS geo_field,
    ST_SRID(geometry) AS srid,
    'shp' AS cte
FROM
    shp
GROUP BY 
1,2,3

This threw the same original error: enter image description here

I'm not sure where to go from here, any help would be most appreiated!

1

There are 1 answers

0
JGH On

The srid in sd_geo_uk_counties is not set, so 0 is assumed and it is an invalid value for st_transform.

Try setting it first:

SELECT 
    ctyua23nm AS county,
    ST_Transform(st_srid(geometry,27700), 4326) AS geometry
FROM 
    sd_geo_uk_counties

Idealy you would set the srid at the column level and update (once) every geometry to set the srid.

On a side note, you have swapped latitude and longitude when building the point, it should be longitude first: ST_SetSRID(ST_Point(longitude, latitude),4326) AS lat_long