# Create polygon in LibGEOS.jl

**URL:** <https://discourse.julialang.org/t/create-polygon-in-libgeos-jl/93483>\
**Category:** General Usage\
**Tags:** polygons, libgeos\
**Created:** [January 25, 2023, 2:24am UTC](https://discourse.julialang.org/t/create-polygon-in-libgeos-jl/93483 "2023-01-25T02:24:57Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![SDA](https://avatars.discourse-cdn.com/v4/letter/s/41988e/32.png) [@SDA](https://discourse.julialang.org/u/SDA)\
**Post date:** [January 25, 2023, 2:24am UTC](https://discourse.julialang.org/t/create-polygon-in-libgeos-jl/93483/1 "2023-01-25T02:24:57Z")

</div>

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

Any clues?

---

<div class="post-metadata">

**Author:** ![joa-quim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joa-quim/32/227_2.png) [@joa-quim](https://discourse.julialang.org/u/joa-quim)\
**Post date:** [January 25, 2023, 2:54am UTC](https://discourse.julialang.org/t/create-polygon-in-libgeos-jl/93483/2 "2023-01-25T02:54:52Z")

</div>

You may want to try GMT that also wraps GDAL and has an IMO easier interface. The [intersection](https://www.generic-mapping-tools.org/GMTjl_doc/documentation/gdalfuns/#intersection) can be run using `Mx2` matrices directly.

```julia
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

```

---

<div class="post-metadata">

**Author:** ![SDA](https://avatars.discourse-cdn.com/v4/letter/s/41988e/32.png) [@SDA](https://discourse.julialang.org/u/SDA)\
**Post date:** [January 31, 2023, 4:51am UTC](https://discourse.julialang.org/t/create-polygon-in-libgeos-jl/93483/3 "2023-01-31T04:51:59Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![joa-quim](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/joa-quim/32/227_2.png) [@joa-quim](https://discourse.julialang.org/u/joa-quim)\
**Post date:** [January 31, 2023, 11:59am UTC](https://discourse.julialang.org/t/create-polygon-in-libgeos-jl/93483/4 "2023-01-31T11:59:07Z")

</div>

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

```julia
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)

```

 ![GMTjl_tmp](https://global.discourse-cdn.com/julialang/original/3X/0/e/0edf99def7b72da038848f756b79a065b556936d.png)

---

<div class="post-metadata">

**Author:** ![PeterSimon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petersimon/32/25193_2.png) [@PeterSimon](https://discourse.julialang.org/u/PeterSimon)\
**Post date:** [January 31, 2023, 7:35pm UTC](https://discourse.julialang.org/t/create-polygon-in-libgeos-jl/93483/5 "2023-01-31T19:35:39Z")

</div>

I did some experimenting with LibGEOS in an IJulia notebook a while ago. I’ve uploaded the notebook to [this gist](https://gist.github.com/simonp0420/256412e3bd211ffea4aadf7a64f041db) where I think you can see some examples of how do what you’ve asked.