0

I am storing OS Open Local data, as downloaded here: https://www.ordnancesurvey.co.uk/opendatadownload/products.html (top result).

I've imported the .shp files.

Running the .prj through here: http://prj2epsg.org/search produces:

    PROJCS["British_National_Grid",GEOGCS["GCS_OSGB_1936",DATUM["D_OSGB_1936",
SPHEROID["Airy_1830",6377563.396,299.3249646]],PRIMEM["Greenwich",0.0],
UNIT["Degree",0.0174532925199433]],PROJECTION["Transverse_Mercator"],
PARAMETER["False_Easting",400000.0],PARAMETER["False_Northing",-100000.0],
PARAMETER["Central_Meridian",-2.0],PARAMETER["Scale_Factor",0.9996012717],
PARAMETER["Latitude_Of_Origin",49.0],UNIT["Meter",1.0],AUTHORITY["EPSG",27700]]

So in my postgis db I run:

SELECT srid FROM spatial_ref_sys WHERE srtext LIKE '%British_National_Grid%';

Which produces two results: 27700 and 7405

By inspection the former corresponds exactly to my .prj

Verifiying my SRID on my column:

SELECT Find_SRID('', 'buildings', 'geometry');

 find_srid 
-----------
         0
(1 row)

Fine, so casting the column to 27700 should do the trick:

SELECT *
FROM buildings AS b
WHERE ST_DWithin(ST_SetSRID(b.geometry, 27700),
ST_Transform(ST_SetSRID(ST_MakePoint(-0.1047, 51.5236), 4326),27700),
500);

(Note that (51.5246, -0.1047) is downtown London)

This query produces no results. interestingly, increasing the radius to 15000, produces 8748 results. Casting an arbitrary one to 4326 and using ST_AsGeoJSON to produce a readable result gives a location in the suburbs (northwest) of London which suggests that my coordinates are somewhat off...

EDIT: Flipping the coordinates to correspond to more traditional coordinate systems also generates no results (more or less at any radius, I've gone up to 500000):

SELECT *
FROM buildings AS b
WHERE ST_DWithin(ST_SetSRID(b.geometry, 27700),
ST_Transform(ST_SetSRID(ST_MakePoint(51.5236, -0.1047), 4326),27700),
500000);

How can I figure out where and/or how things are going wrong? And how can then fix this?

  • Question edited to include result of flipping coordinate system to be correct – Abraham P Jun 16 '17 at 10:41
  • you need to set the srid of your buildings layer to 27700 not 0 – Ian Turton Jun 16 '17 at 10:42
  • isn't that being done by the ST_SetSRID(b.geometry, 27700)? – Abraham P Jun 16 '17 at 10:49
  • {51.5246, -0.1047} is in the middle of the Atlantic Ocean, near the equator. – Vince Jun 16 '17 at 10:52
  • 1
    Coordinates are entered in X,Y order, so longitude,latitude is proper. – Vince Jun 16 '17 at 11:03
2

I think you may not have imported your shapefile correctly as this works for me.

First create the table (as shp2psql had an issue with the add geom column):

CREATE TABLE "ian"."buildings" (gid serial,
"id" varchar(38),
"featcode" int4);
ALTER TABLE "ian"."buildings" ADD PRIMARY KEY (gid);
SELECT AddGeometryColumn('ian','buildings','geom','27700','MULTIPOLYGON',2,true);

then import the shapefile:

shp2pgsql -s 27700 -I -a Building.shp ian.buildings | psql -f -

Then the query

SELECT *
FROM buildings AS b
WHERE ST_DWithin(ST_SetSRID(b.geom, 27700),
ST_Transform(ST_SetSRID(ST_MakePoint(-0.1047,51.5236), 4326),27700),
5000);

returns 7180 buildings.

Note I was wrong in my assumption that PostGIS would do the right thing with lat/lon v. lon/lat as can be seen with the following test.

SELECT ST_ASTEXT(ST_Transform(ST_SetSRID(ST_MakePoint(-0.1047,51.5236), 4326),27700));

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.