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?
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 :
ST_Centroid
to get the centroid of the geometry - a pointST_Y
function to obtain the latitudeST_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!