Tuesday, 9 April 2019

postgis - How can I split a polygon into two equal parts along a N-S axis?


I have a PostGIS database with a table called 'sites' that includes a field sitename and a field geometry.



If I select, e.g.


select sitename, ST_AsGeoJSON(geometry) from sites limit 1;

I get:



Plot 7 | {"type":"Polygon","coordinates": [[[-111.974951718315,33.0745275800293,353], [-111.974935360975,33.0745276003588,353], [-111.974935370263,33.0745546649121,353], [-111.974951727608,33.0745546445826,353],[-111.974951718315,33.0745275800293,353]]]}



These are rectangles oriented very close to N-S. I would like to divide each plot in half along (or close to) its N-S axis. Is there an easy way to do this?


Solutions in PostGIS-SQL preferred. Ultimately I plan to write insert statements like


insert into sites (name, geometry) 

values ('Site 1 E', (select ??east half of geometry?? from sites where name = 'Site 1'))

The closest I've gotten is ST_Split, which requires the polygon as well as the 'blade' (line used to split), but it is not clear how I can compute the blade from the geometry.



Answer




I would like to divide each plot in half along (or close to) its N-S axis. Is there an easy way to do this? [...] can compute the blade from the geometry.



Something like?


xmin+((xmax-xmin)/2)


Here we create a sample data set foo. In it we create a polygon, we calculate the xaxis to bisect the polygon with, and we create the blade.


Then we do the bisection and display the results with ST_AsText


WITH foo AS (
SELECT geom, blade
FROM ST_MakeEnvelope(-10,-10,10,10) AS geom
CROSS JOIN LATERAL ( SELECT ST_xMin(geom) + (ST_xMax(geom) - ST_xMin(geom)) / 2 ) AS axis(x)
CROSS JOIN LATERAL ST_MakeLine(
ST_MakePoint(axis.x, ST_yMin(geom)),
ST_MakePoint(axis.x, ST_yMax(geom))
) AS blade

)
SELECT ST_AsText(geom),
ST_AsText(blade),
ST_AsText( ST_Split(geom, blade) )
FROM foo;
st_astext | st_astext | st_astext
------------------------------------------------+------------------------+----------------------------------------------------------------------------------------------------------
POLYGON((-10 -10,-10 10,10 10,10 -10,-10 -10)) | LINESTRING(0 -10,0 10) | GEOMETRYCOLLECTION(POLYGON((-10 -10,-10 10,0 10,0 -10,-10 -10)),POLYGON((0 10,10 10,10 -10,0 -10,0 10)))

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