Create regular point grids in polygons

Q&A GMT.jl Q&A
1

Hi,

I would like to generate a regular grid of points inside a list of polygons.

I did this in Python:

def generate_grid_in_polygon(resolution, polygon):
    points = []
    minx, miny, maxx, maxy = polygon.bounds
    for x in np.arange(minx, maxx, resolution):
        for y in np.arange(miny, maxy, resolution):
            pnt = Point(x, y)
            if polygon.contains(pnt):
                points.append(pnt)
    return points

Is there any way that I can do this using GMT.jl and / or standard Julia Array?

Thanks!

2

One way (if I understood the Python code)

X,Y = meshgrid(1:50, 1:50);

# A triangle polgon
polyg = [10 10; 25 25; 40 10; 10 10];

# Pick only the points inside the triangle
D = gmtselect([X[:] Y[:]], polygon=polyg);

# Show it
plot(D, marker=:point, show=1)

GMTjl_tmp
GMTjl_tmp1756Ɨ1211 26.8 KB

1 Like
3

Thanks - Is there a way to have the bounds of a polygon using GMT.jl to create the meshgrid within the polygon limits before gmtselect?

something like this maybe:

X,Y = meshgrid(minimum(polygon.xcoordinates):maximum(polygon.xcoordinates), 
                          minimum(polygon.ycoordinates):maximum(polygon.ycoordinates) )
4

If you convert a matrix into a GMTdataset with D = mat2ds(mat) the D type has fields D.bbox and D.ds_bbox (when D is a multi-segment a vector of Dā€™s) that holds that info.

GMTdatasets are quite complete types. You should read the GMT types section of the docs

1 Like
5

Thank you very much for all the help - I gain a LOT of time with your advices