2

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

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: Modelbuilder function

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.

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).

  • 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 4 '18 at 7: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. – ThingumaBob Dec 4 '18 at 8: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 5 '18 at 9: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)); – ThingumaBob Dec 5 '18 at 11:39
  • 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; – ThingumaBob Dec 7 '18 at 9:21

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.