I am trying to create a polygon using LibGEOS.jl. The docs give an example for the case where the polygon vertices are input manually:

using LibGEOS
p1 = readgeom(“POLYGON((0 0,1 0,1 1,0 0))”)

I would like to instead create this polygon by referencing a vector that contains the vertices.

For example, I can do this to create a point:

x=1
y=1
readgeom(“POINT($x $y)”) #this works

But I can’t figure out how to extend it to the case for a polygon. Here’s one example I’ve tried that doesn’t work:

xs = 0, 1, 2, 0 #defining this as a vector also doesn’t work
ys = 0, 0, 1, 0
readgeom(“POLYGON($xs $ys)”)

Edited to include: I’ve also tinkered with ArchGDAL though I think this might be overcomplicating things.

using ArchGDAL
xs = [0, 1, 2, 0]
ys = [0, 0, 1, 0]
ArchGDAL.createpolygon(xs,ys)
But this returns Geometry: POLYGON ((0 0, ...)
while I need only the POLYGON((0 0, ...) part (to feed in to LibGEOS.intersection(...) )

Thank you @joa-quim, I wasn’t aware of GMT so this was interesting. It doesn’t quite achieve what I want as it returns only the vertices of the intersection of two shapes. It is fine for most cases but I need to intersect some shapes with a circle and return the entire circle segment or at least an approximation of it.

i.e., I am seeking the intersection of the shapes in the sense of the area common to both shapes, rather than the intersecting points of them.

I did some experimenting with LibGEOS in an IJulia notebook a while ago. I’ve uploaded the notebook to this gist where I think you can see some examples of how do what you’ve asked.