postgresqlamazon-web-servicesamazon-rdspostgisspatial-query

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


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

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:

enter image description here

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?


Solution

  • 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:

    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!

    enter image description here