Home > Net >  How to find parcels that are bordered with my parcel on Postgresql
How to find parcels that are bordered with my parcel on Postgresql

Time:12-12

My problem is: I have 1000 parcels and 1 of them is mine, what I have to do is to find which are my neighbourhoods which means whose parcels are touched with my parcel. The perfect example is showed on the picture that I have sent you here. Red rectangle is my parcel and black rectangles are my neighbourhoods what I have to do is to find who are they. I have created 2 tables, first table(parcels) contains all parcels without mine and second table(Myparcel) contains just my parcel. As spatial data, I have the geometry of all parcels ( column geom) and the number of parcels.My parcel is red color and neighbourhood parcels are black

I'm trying to solve a problem, I have been trying for several days to solve it but no results! So the final results is that I wanna know which parcels are touched with my parcel. I have created 2 tables and here is my code which is not working

select st_touches(k.geom, s.geom) from kultura5_maxsip1 k, kultura5_maxsip s

CodePudding user response:

ST_Touches() likely won't return much, if anything, because typically geometries are naturally shifted around, unless you ST_Snap() them together, or apply a uniform ST_SnapToGrid(). Also, the intuitive understanding of "touch" is misleading here, because ST_Touches() returns True if one of your shapes ends precisely where the other one begins, without crossing or overlapping. If one's entirely inside the other, they aren't considered touching either.

It's good to read about DE-9IM or simply go through the PostGIS doc to see which functions look for what types of spatial relationships. You can easily test all of them constructing very simple shapes with ST_GeomFromText() and then inspecting them with ST_AsText() or in QGIS. Demo:

create table my_parcel(
  id serial primary key, 
  geom geometry(POLYGON,--type of shape
                21781   --SRID, Switzerland in this case, coordinates in METERS
               )
  );
--always build an index on the geometry column to speed things up
create index on my_parcel using gist(geom);

insert into my_parcel(geom) values 
(ST_GeomFromText('POLYGON((0 0,1 0,1 1,0 1,0 0))')); --a 1x1 square, north-east of origin

You can read about SRIDs on EPSG.io. It's important to remember that they define your target area limits but also the unit of measurement. If you use global 4326 and then measure an ST_Distance() between shapes in that system, it'll be in degrees, not meters. If you use 21781, you won't be able to encode coordinates outside Switzerland with some buffer. There are systems in miles/feet, also.

create table all_parcels(like my_parcel including all);

insert into all_parcels(geom) 
select st_translate(geom,move.x,move.y) --generates a copy moved by an offset
from my_parcel, 
  (values --these offsets define chess king movements, 1 in all directions
  (0,1),  --move up             /north 
  (1,1),  --move right and up   /north-east
  (1,0),  --move right          /east
  (-1,0), --move left           /west
  (-1,-1),--move left down      /south-west
  (0,-1), --move down           /south
  (1,-1)  --move right and down /south-east
  ) as move(x,y);

These parcels will touch in ST_Touches() sense, because they re-use the same, perfectly round coordinates.

select id,st_astext(geom) from all_parcels;
-- id |              st_astext
------ --------------------------------------
--  2 | POLYGON((0 1,1 1,1 2,0 2,0 1))
--  3 | POLYGON((1 1,2 1,2 2,1 2,1 1))
--  4 | POLYGON((1 0,2 0,2 1,1 1,1 0))
--  5 | POLYGON((-1 0,0 0,0 1,-1 1,-1 0))
--  6 | POLYGON((-1 -1,0 -1,0 0,-1 0,-1 -1))
--  7 | POLYGON((0 -1,1 -1,1 0,0 0,0 -1))
--  8 | POLYGON((1 -1,2 -1,2 0,1 0,1 -1))
select b.id
from my_parcel a,
  all_parcels b
where st_touches(a.geom,b.geom);
-- id
------
--  2
--  3
--  4
--  5
--  6
--  7
--  8

Real-life data usually comes from error-prone measurements in suboptimal conditions. As soon as they move a tiniest bit, they no longer touch:

create table all_parcels_slightly_shifted(like all_parcels including all);

insert into all_parcels_slightly_shifted 
select 
  id, 
  st_translate(geom, random()/100, random()/100) --random milimeter offset
from all_parcels;
-- id |                                                                                             st_astext                                            
------ ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
--  2 | POLYGON((0.005739247141314 1.002288654993983,1.005739247141314 1.002288654993983,1.005739247141314 2.002288654993983,0.005739247141314 2.002288654993983,0.005739247141314 1.002288654993983))
--  3 | POLYGON((1.005088105236068 1.000882987726616,2.005088105236068 1.000882987726616,2.005088105236068 2.000882987726616,1.005088105236068 2.000882987726616,1.005088105236068 1.000882987726616))
--  4 | POLYGON((1.001129227264427 0.007497058517941,2.001129227264427 0.007497058517941,2.001129227264427 1.007497058517941,1.001129227264427 1.007497058517941,1.001129227264427 0.007497058517941))
--  5 | POLYGON((-0.991909588170889 0.006005249136691,0.008090411829111 0.006005249136691,0.008090411829111 1.006005249136691,-0.991909588170889 1.006005249136691,-0.991909588170889 0.006005249136691))
--  6 | POLYGON((-0.993659995690557 -0.99898573996986,0.006340004309443 -0.99898573996986,0.006340004309443 0.00101426003014,-0.993659995690557 0.00101426003014,-0.993659995690557 -0.99898573996986))
--  7 | POLYGON((0.009123227228849 -0.992001817630975,1.009123227228849 -0.992001817630975,1.009123227228849 0.007998182369025,0.009123227228849 0.007998182369025,0.009123227228849 -0.992001817630975))
--  8 | POLYGON((1.005054414974922 -0.996745585497164,2.005054414974923 -0.996745585497164,2.005054414974923 0.003254414502835,1.005054414974922 0.003254414502835,1.005054414974922 -0.996745585497164))
select b.id 
from my_parcel a, 
  all_parcels_slightly_shifted b 
where st_touches(a.geom,b.geom);
--<nothing>

But if you ask which ones now intersect (have any spatial relationship) or are within a certain distance (that you'd consider close enough to agree they are touching in the colloquial sense):

select b.id, st_astext(b.geom)
from my_parcel a,
  all_parcels_slightly_shifted b
where st_intersects(a.geom,b.geom);
-- id |                                                                                             st_astext
------ ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
--  5 | POLYGON((-0.991909588170889 0.006005249136691,0.008090411829111 0.006005249136691,0.008090411829111 1.006005249136691,-0.991909588170889 1.006005249136691,-0.991909588170889 0.006005249136691))
--  6 | POLYGON((-0.993659995690557 -0.99898573996986,0.006340004309443 -0.99898573996986,0.006340004309443 0.00101426003014,-0.993659995690557 0.00101426003014,-0.993659995690557 -0.99898573996986))
--  7 | POLYGON((0.009123227228849 -0.992001817630975,1.009123227228849 -0.992001817630975,1.009123227228849 0.007998182369025,0.009123227228849 0.007998182369025,0.009123227228849 -0.992001817630975))
select b.id, st_astext(b.geom)
from my_parcel a,
  all_parcels_slightly_shifted b
where st_dwithin(a.geom,b.geom,0.01);
-- id |                                                                                             st_astext
------ ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
--  2 | POLYGON((0.005739247141314 1.002288654993983,1.005739247141314 1.002288654993983,1.005739247141314 2.002288654993983,0.005739247141314 2.002288654993983,0.005739247141314 1.002288654993983))
--  3 | POLYGON((1.005088105236068 1.000882987726616,2.005088105236068 1.000882987726616,2.005088105236068 2.000882987726616,1.005088105236068 2.000882987726616,1.005088105236068 1.000882987726616))
--  4 | POLYGON((1.001129227264427 0.007497058517941,2.001129227264427 0.007497058517941,2.001129227264427 1.007497058517941,1.001129227264427 1.007497058517941,1.001129227264427 0.007497058517941))
--  5 | POLYGON((-0.991909588170889 0.006005249136691,0.008090411829111 0.006005249136691,0.008090411829111 1.006005249136691,-0.991909588170889 1.006005249136691,-0.991909588170889 0.006005249136691))
--  6 | POLYGON((-0.993659995690557 -0.99898573996986,0.006340004309443 -0.99898573996986,0.006340004309443 0.00101426003014,-0.993659995690557 0.00101426003014,-0.993659995690557 -0.99898573996986))
--  7 | POLYGON((0.009123227228849 -0.992001817630975,1.009123227228849 -0.992001817630975,1.009123227228849 0.007998182369025,0.009123227228849 0.007998182369025,0.009123227228849 -0.992001817630975))
--  8 | POLYGON((1.005054414974922 -0.996745585497164,2.005054414974923 -0.996745585497164,2.005054414974923 0.003254414502835,1.005054414974922 0.003254414502835,1.005054414974922 -0.996745585497164))

You can read about Hausdorff distance, which is a bit better at telling you how close entire shapes are to each other, not just how long is the shortest shortcut between them. ST_DWithin() is True based on whether the minimum distance, the shortest shortcut between shapes is below the limit.

CodePudding user response:

You could try using DBScan (implemented in PostGIS) for this use case to determine neighbouring clusters. Please give this a shot:

SELECT *, 
       ST_ClusterDBSCAN(geom, 0, 1) OVER() AS cluster_id 
FROM <your_table>;

cluster_id represents the cluster a row belongs to. With a epsilon value of 0, the function clusters by intersection only. If your polygons are not intersecting, you could buffer them by a small value.

  • Related