Create polygon in LibGEOS.jl

General Usage
,
1

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?

2

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
3

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.

4

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)

GMTjl_tmp
GMTjl_tmp1800×1508 43.9 KB

1 Like
5

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