Tuesday, 17 May 2016

srid - Control number of segments in ST_Buffer for geographies?


I've been attempting to construct circles centered around a list of points in PostGIS with good success. I've opted to circumvent having to look up UTM zones and using optimal SRID by simple letting ST_Buffer do the work for me by supplying a geography rather than a geometry.


However, in doing so I lose the ability to choose the number of segments for the resulting circle geometry, which defaults to 33 (not sure how this number is calculated, I make circle areas of radius 20 km). Is there any way to circumvent this and choose the number of segments myself? For my query of points I use:


update 
set poly = geometry(st_buffer(geography(point), 20000));


The output geometry will default to SRID 4326.


 
Answer
 

Buffer does not seem to support any extra parameters with geography http://postgis.net/docs/ST_Buffer.html. However, with geometry you can control the number of segments with num_seg_quarter_circle.

You need to cast your geography into geometry before using ST_Buffer with num_seg_quarter_circle parameter and finally cast the result back to geography. Try this query for creating a geography which is an input point buffered by 1000 units with 10 segments per quadrant.

SELECT ST_Buffer(
ST_GeographyFromText(
'SRID=4326;POINT(-110 30)')::geometry,1000,10)::geography;

This query works but the result is nonsense because cast from geography into geometry yields "SRID=4326;POINT(-110 30)". 1000 units buffer for EPSG:4326 means a buffer of 1000 degrees. The query must be enhanced to include ST_Transform for projecting the data into some reasonable projected coordinate system using meters as unit before applying ST_Buffer. After buffering reverse ST_Transform is needed.


A prototype query for a point at (10,10) lon-lat going via EPSG:32630 which is WGS84 UTM zone 30:

select ST_Transform(
ST_Buffer(
ST_Transform(ST_GeographyFromText('SRID=4326;POINT(10 10)')::geometry,32630),1000,10),4326)::geography

Note: The second ST_Transform from EPSG:32630 into EPSG:4326 is probably not necessary because cast to geography takes care of that automatically. However, I feel that it is easier to understand the process with both forward and reverse ST_Transform.







-



































No comments:















Post a Comment

















Subscribe to:
Post Comments (Atom)







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

Blog Archive