Select layer by location in PostgreSQL

Question

I have a footprint of some buildings (see the picture).

The selected buildings have an area below 30 sqm. As you can in the picture see some buildings share almost borders with a building bigger than 30 sqm. When others don't. The building below 30 sqm but shares a border with a bigger building I want to include in my calculation. The others below 30 sqm should be sorted away.

I have succeeded with this through desktop GIS, and the solution looks like this:

But I want to do the same process in postgres.

I can select the buildings pending on their area and I have put them in to two new tables:

Create TABLE building_below as Select * from building where area<30 and Create TABLE building_above as Select * from building where area>30

I can not figure out how to select layer by location in PostgreSQL.

I have tried with inspiration from: Join polygons by location PostgreSQL/PostGIS Merging adjacent polygons in shapefile that has been split at tile boundaries?

But it didn't work.

I have tried this:

SELECT DISTINCT ON (s.gid) s.gid,*

 FROM building_below s

LEFT JOIN building_above h ON ST_DWithin(s.geom, h.geom, 0.01)

ORDER BY s.gid, ST_Distance(s.geom, h.geom);

But it just select all the rows/buildings in building_below.

I know from desktop GIS, that it should not do that.

geozelot · Accepted Answer · 2018-12-07T09:26:30.193

2

Try

SELECT a.*
FROM  building AS a
JOIN  building AS b
  ON  ST_Intersects(a.geom, b.geom)
WHERE a.area <= 30
  AND b.area > 30;

Make sure you have a spatial index in place and run VACUUM ANALYZE building; prior to execution!

If your buildings only almost share a border, i.e. don't actually intersect, use

... ST_DWithin(a.geom, b.geom, <threshold>) ...

instead.

Note that the latter will need your geom either projected or in type GEOGRAPHY to make sense (i.e. if your geom is in EPSG:4326 or any geographic CRS, units will be in degrees and thus rather useless. You can then, if in EPSG:4326, however, easily cast to GEOGRAPHY to implicitly use meter as unit).

edited Dec 07 '18 at 09:26

answered Dec 03 '18 at 15:01

geozelot

30,050
4
32
56

My Geom is in 4326. But the area is calculated in 25832. What should do then? When I run this:
SELECT * FROM bygning AS B JOIN bygning AS D ON public.ST_Dwithin(B.geom, D.geom, 0.5) WHERE B.area_25832 <= 30 AND D.area_25832 > 30;

It returns +300 rows, and their is only 60 polygons.
– N.Juul Dec 04 '18 at 07:56
ST_DWithin(B.geom::GEOGRAPHY, D.geom::GEOGRAPHY, 0.5) ... will find geometries within 0.5 meter. problem now is that the index on GEOMETRY cannot be used; either create a functional spatial index ... ON (CAST(geom, GEOGRAPHY)) ... on both tables, create one ... ON (ST_Transform(geom, 25832)) ... on both tables and transform on-the-fly, or reproject (or retype) and reindex both tables into EPSG:25832 (or GEOGRAPHY). you could try without all that and using 0.000005 (degrees) maybe...but that's a last resort. – geozelot Dec 04 '18 at 08:28
I have tried the following sentence; <--alter table bygning add column geom_25832 geography; update bygning set geom_25832=(ST_Transform(geom, 25832))>
the output is: ERROR: Only lon/lat coordinate systems are supported in geography. I guess, i am not sure, what you want me to do? – N.Juul Dec 05 '18 at 09:22
@N.Juul the GEOGRAPHY type needs the geometries in EPSG:4326 (-> WGS 84 lat/lons): if your initial geometries are EPSG:4326, I'd suggest to just add a functional index with a cast: CREATE INDEX bygning_geog_idx ON bygning USING GIST (CAST(geom, GEOGRAPHY)); – geozelot Dec 05 '18 at 11:39
I have tried this:
drop index bygning_geog_idx; CREATE INDEX bygning_geog_idx ON bygning USING GIST (CAST(geom as GEOGRAPHY));

SELECT * FROM bygning AS a JOIN bygning AS b ON ST_Dwithin(a.geom, b.geom) WHERE a.area <= 30 AND b.area > 30;

And I have run the vacuum analyze in prior.

But it still does not work.

I have one table buildings, with a geom in EPSG:4326. I have two tables with buidlings below 30 sqm and above 30 sqm.

Buildings below 30 sqm, but within 0,5 meters from a building above 30 sqm, should also be a part of the calculation.

How to do that?
– N.Juul Dec 06 '18 at 08:42
@N.Juul assuming your geom column is still type GEOMETRY and referenced in EPSG:4326, and you did add that index, you'd have to then cast the geom to GEOGRAPHY: --> SELECT * FROM bygning AS a JOIN bygning AS b ON ST_DWithin(a.geom::GEOGRAPHY, b.geom::GEOGRAPHY, 0.5) WHERE a.area <= 30 AND b.area > 30; – geozelot Dec 06 '18 at 17:03

With your help, i have make it work - almost! I select the buildings thats below 30 sqm, and within 0,5m from a building above 30 sqm. I want to put this into a table with buildings above 30sqm - this should be pretty easy, but it does not work.

  'I have this code:
  --create table bygning_30 as
  --select * from bygning where area_25832>30;
  insert into bygning_30
  SELECT *
  FROM  bygning AS a
 JOIN  bygning AS b
 ON  ST_DWithin(a.geom::GEOGRAPHY, b.geom::GEOGRAPHY, 0.5)
 WHERE a.area_25832 <= 30 
 and b.area_25832 >30

– N.Juul Dec 07 '18 at 09:15

1

@N.Juul what's the problem? any errors? how about CREATE TABLE bygning_30 AS SELECT * FROM bygning WHERE area_25832 > 30 UNION ALL SELECT a.* FROM bygning AS a JOIN bygning AS b ON ST_DWithin(a.geom::GEOGRAPHY, b.geom::GEOGRAPHY, 0.5) WHERE a.area_25832 <= 30 AND b.area_25832 > 30; – geozelot Dec 07 '18 at 09:21
In my solution it says the error is: ERROR: INSERT has more expressions than target columns LINE 1: insert into bygning_30 SELECT * FROM bygning AS a JOIN bygni... When I run your solution the error is:
ERROR: each UNION query must have the same number of columns LINE 4: UNION ALL SELECT * FROM bygning AS a JOIN bygning AS b
– N.Juul Dec 07 '18 at 09:25
@N.Juul ...ah yes, I forgot one thing: you will need to specify the table from which to select the columns! in your example (and mine thus far), SELECT * FROM ... will select all columns from both tables of the join, i.e. it duplicates the same columns in this case of a self-join. I updated the answer, and the query in the comment I just wrote, with SELECT a.* FROM .... – geozelot Dec 07 '18 at 09:29
Perfect! It works - thanks for your patience and help!!! My project group appreciates that! – N.Juul Dec 07 '18 at 09:35

Select layer by location in PostgreSQL

1 Answers1