Why are UK county geometries showing up in the Gulf of Guinea when using PostGIS?

Question

I have two tables in a AWS Postgres server (I am querying using DBeaver).

• Table 1 is a shapefile which I imported from an S3 bucket called uk_counties, which lists a geometry for all counties in the UK (SRID = 27700 - shapefile can be found here).

• Table 2 is called demand_origin and contains the fields origin_city (places in the UK) and two fields latitude and longitude (which contains the latitutde and logitudes of places in origin_city).

By joining uk_counties (the shapefile) and demand_origin, I want to assign a county to the places in the UK by finding where their respective lat longs intersect with the shapefiles for the counties.

The problem is that this join kept appearing empty. After much investigation, it appears that the geometries in the shapefile are mapped to the Gulf of Guinea (in and around lat = 0, long = 0).

See shapefile for Leicester as an example:

My guess here is that the geometries are using northings/westings as opposed to lat/longs, hence PostGIS is mapping these shapes to around lat=0, long=0, and as such potentially needs to be recalibrated.

Here is a snippet of what the geometries look like for Leicester and other counties:

ctyua23nm	geometry
York	POLYGON ((464216.99650000036 462110.9015999995, 464266.20299999975 462093.0965999998, 464270.3008000003 462094.3962999992, 464289.8992999997 4
Derby	POLYGON ((434972.3005999997 341311.7999000009, 434986.3008000003 341310.40029999986, 435000.0536000002 341314.6991000008, 435015.29860000033 3
Cambridgeshire	MULTIPOLYGON (((529832.0997000001 300053.7999000009, 529880.2013999997 300039.1041000001, 529903.7987000002 300036.6048000008, 529914.20399999

What's the issue?

Answer 1

FWIW, when I downloaded the Shapefile from the link above, and imported it like so:

shp2pgsql -I -s 27700 -g geometry CTYUA_MAY_2023_UK_BFC.shp uk_counties | psql -d mydb

I obtain a much different schema:

\d+ uk_counties

   Column   |             Type             | Collation | Nullable |                 Default                  | Storage  | Compression | Stats target | Description
------------+------------------------------+-----------+----------+------------------------------------------+----------+-------------+--------------+-------------
 gid        | integer                      |           | not null | nextval('uk_counties_gid_seq'::regclass) | plain    |             |              |
 ctyua23cd  | character varying(9)         |           |          |                                          | extended |             |              |
 ctyua23nm  | character varying(36)        |           |          |                                          | extended |             |              |
 ctyua23nmw | character varying(24)        |           |          |                                          | extended |             |              |
 bng_e      | double precision             |           |          |                                          | plain    |             |              |
 bng_n      | double precision             |           |          |                                          | plain    |             |              |
 long       | double precision             |           |          |                                          | plain    |             |              |
 lat        | double precision             |           |          |                                          | plain    |             |              |
 globalid   | character varying(38)        |           |          |                                          | extended |             |              |
 geometry   | geometry(MultiPolygon,27700) |           |          |                                          | main     |             |              |
Indexes:
    "uk_counties_pkey" PRIMARY KEY, btree (gid)
    "uk_counties_geometry_idx" gist (geometry)
Access method: heap

And I can easily get a correct latitude & longitude value:

SELECT
    ctyua23nm,
    long,
    lat
FROM
    uk_counties
WHERE
    ctyua23nm = 'Leicester';

 ctyua23nm |  long   |   lat
-----------+---------+---------
 Leicester | -1.1304 | 52.6359
(1 row)

It might easier for you to just reimport the Shapefile, as I'm not sure how you've imported it.

If you're stuck with your current schema, please continue.

---

the geometries are using northings/westings as opposed to lat/longs

Correct.

SRID 27700 / EPSG:27700 aligns to the Ordnance Survey National Grid (OSGB) or more commonly known as the British National Grid (BNG). It basically uses projected eastings and northings, which are in metres - not latitude and longitude, which are in degrees.

You can see this by running:

SELECT srtext FROM spatial_ref_sys WHERE srid = 27700;

PROJCS["OSGB 1936 / British National Grid",
  GEOGCS["OSGB 1936",
    DATUM["OSGB_1936",
      SPHEROID["Airy 1830", 6377563.396, 299.3249646,
        AUTHORITY["EPSG", "7001"]],
      TOWGS84[446.448, -125.157, 542.06, 0.15, 0.247, 0.842, -20.489],
      AUTHORITY["EPSG", "6277"]],
    PRIMEM["Greenwich", 0, AUTHORITY["EPSG", "8901"]],
    UNIT["degree", 0.0174532925199433, AUTHORITY["EPSG", "9122"]],
    AUTHORITY["EPSG", "4277"]],
  PROJECTION["Transverse_Mercator"],
  PARAMETER["latitude_of_origin", 49],
  PARAMETER["central_meridian", -2],
  PARAMETER["scale_factor", 0.9996012717],
  PARAMETER["false_easting", 400000],
  PARAMETER["false_northing", -100000],
  UNIT["metre", 1, AUTHORITY["EPSG", "9001"]],
  AXIS["Easting", EAST],
  AXIS["Northing", NORTH],
  AUTHORITY["EPSG", "27700"]]

Note UNIT["metre", AXIS["Easting", EAST] and AXIS["Northing", NORTH].

You most likely want to project your data to use SRID 4326 / EPSG:4326, which aligns to WGS84 i.e. the standard for GPS.

It's a 2D geographic coordinate system i.e. uses latitude and longitude.

---

To project your SRID 27700 table to latitude and longitude coordinates (SRID 4326), firstly verify that the SRID for your geometry column is set to SRID 27700 using Find_SRID.

SELECT
    Find_SRID ('public', 'uk_counties', 'geometry');

The output should be:

 find_srid
-----------
     27700
(1 row)

If not, set the SRID for the column using UpdateGeometrySRID:

SELECT
    UpdateGeometrySRID ('uk_counties', 'geometry', 27700);

                updategeometrysrid
---------------------------------------------------
 public.uk_counties.geometry SRID changed to 27700
(1 row)

Then, use ST_Transform function in PostGIS to do the conversion:

Returns a new geometry with its coordinates transformed to a different spatial reference system.

ALTER TABLE uk_counties
    ALTER COLUMN geometry TYPE geometry(geometry, 4326)
    USING ST_Transform (geometry, 4326);

ALTER TABLE

Once done, we can then get latitude & longitudes from the geometry by :

1. using ST_Centroid to get the centroid of the geometry - a point
2. apply the ST_Y function to obtain the latitude
3. apply the ST_X function to obtain the longitude.

Let's create 2 new columns to store the data:

ALTER TABLE uk_counties
    ADD COLUMN latitude DOUBLE PRECISION,
    ADD COLUMN longitude DOUBLE PRECISION;

ALTER TABLE

And then, filling the columns...

UPDATE
    uk_counties
SET
    latitude = ST_Y (ST_Centroid (geometry)),
    longitude = ST_X (ST_Centroid (geometry));

UPDATE 218

Now, you have 2 new columns:

• latitude
• longitude

The new columns should have correct values:

SELECT
    ctyua23nm,
    latitude,
    longitude
FROM
    uk_counties
LIMIT 10;

          ctyua23nm          |     latitude       |     longitude
-----------------------------+--------------------+---------------------
 Hartlepool                  | 54.669449373886266 | -1.2594048655284926
 Middlesbrough               |  54.54200444097641 | -1.2222619068734741
 Redcar and Cleveland        |  54.55165165356056 | -1.0206837926504952
 Stockton-on-Tees            |  54.56161339847885 | -1.3322804311895529
 Darlington                  | 54.548705061498694 | -1.5526287937957188
 Halton                      | 53.345131982125686 |  -2.712197804308062
 Warrington                  | 53.399109236880015 | -2.5627669775597606
 Blackburn with Darwen       | 53.692643404310815 |  -2.467582012742445
 Blackpool                   |  53.81779799751382 |  -3.029241917065863
 Kingston upon Hull, City of |  53.76405945136756 |  -0.334854807051097
(10 rows)

---

And last but not least, to come full circle with the Leicester example...

SELECT
    ctyua23nm,
    latitudea,
    longitudea
FROM
    uk_counties
WHERE
    ctyua23nm = 'Leicester';

 ctyua23nm |     latitude      |     longitude
-----------+-------------------+---------------------
 Leicester | 52.63792039053044 | -1.1284845708542661
(1 row)

Here's Leicester!