# Create polygon in LibGEOS.jl

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

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

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(...)` )

Any clues?

You may want to try GMT that also wraps GDAL and has an IMO easier interface. The intersection can be run using `Mx2` matrices directly.

``````xs = [0, 1, 2, 0];
ys = [0, 0, 1, 0];

# The intersection with itself returns ... itself
intersection([xs ys], [xs ys])
BoundingBox: [0.0, 2.0, 0.0, 1.0]
4×2 GMTdataset{Float64, 2}
Row │   col.1    col.2
│ Float64  Float64
─────┼──────────────────
1 │     0.0      0.0
2 │     2.0      1.0
3 │     1.0      0.0
4 │     0.0      0.0
``````

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.

Not sure 100% sure I understood what you want, but maybe this:

``````using GMT

circ = arccircle(1, 0, 1, 0, 360);
square = [-0.5 -0.5; -0.5 0.5; 0.5 0.5; 0.5 -0.5; -0.5 -0.5]
arc = gmtselect(circ, polygon=square);
plot(circ, region=(-1.0,2.5,-1.5,1.5), aspect=:equal)
plot!(square)
plot!(arc, lc=:red, show=true)
``````

1 Like

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.

1 Like