Thursday 3 January 2019

postgresql - Adding vertices to polygon in PostGIS?


For example I have rectangular polygon (4 vertices). I want to add mid-point on each line of the polygon, so in total now I have 8 vertices. And from them I want to build new polygon who will have same exterior footprint, but have 8 vertices: enter image description here



What I did and it's not working (I want to do this on every polygon in my buildings.shp): (pseudo algorithm, I'm working with "WITH... AS")




  1. Get all original vertices (The 4 black) ST_DumpPoints on buildings twice(!)


    Prepare borders of the multi-polygon buildings as line-strings in order to apply ST_Line_Interpolate_Point-




  2. All possible lines between same polygon vertices: ST_MakeLine(points1.geom, points2.geom) WHERE (points1.path <> points2.path) AND (points1.id = points2.id)



  3. Buildings borders: ST_ExteriorRing(the_geom) ON (ST_Dump(build.geom)).geom AS the_geom



  4. From lines @ 2 select only the border lines: St_touches(build.geom, lines.geom) AND St_touches(border.geom, lines.geom)=false AND lines.id = build.id AND lines.id = border.id --all points from same building, all lines from same building have same id belonging to the building




  5. Get the mid points: ST_Line_Interpolate_Point(border_lines.geom, 0.5) (green points)




  6. UNION points @ 1 and points @ 5





Now the problem is that I don't have any order in the points @ 6 so it's almost impossible to use ST_MakePolygon, and tricks with ST_ConcaveHull & ST_ConvexHull also not working because they ignores the mid point and eventually I return to where I started.



Answer



Use ST_Snap http://postgis.net/docs/ST_Snap.html as the final step. As input geometries you need your original polygon and extra vertices to be added as MultiPoint. The third term of the function is snapping tolerance.


enter image description here


SELECT ST_Snap(
ST_GeomFromText('POLYGON (( 7 7, 7 11, 11 11, 11 7, 7 7 ))'),
ST_GeomFromtext('MULTIPOINT (( 11 9 ), ( 7 9 ), ( 9 11 ), ( 9 7 ))'),0.01);

Result is a polygon that has new vertices added from the MultiPoint:




POLYGON (( 7 7, 7 9, 7 11, 9 11, 11 11, 11 9, 11 7, 9 7, 7 7 ))



enter image description here


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...