How to check if my gridded X and Y are inside a polygon?

And how to you generate one of those using ndgrid? I’m not saying you never need to store triple coordinates in triple matrices, I’m saying you (should) never need to store a vector as a grid matrix.

Would like to make a plot like the one I made above (forget the hole) with your recipe?
(GMT does NOT gridify the coordinates, and that since ~30 years ago)

I’m not sure I understand the question. It seems like GMT has a pretty good recipe for plotting geographical data, so I don’t see why you would want to reinvent the wheel? Especially if it doesn’t require gridification, there is really no benefit (also I don’t feel like installing 4 GB for a visualization exercise).

Actually, plotting is not my main goal here. I need to produce some NetCDF files with the gridded data so that others can feed the data into their models, for example.

Anyway, your solution solved my problem nicely! :+1:

You could also reshape the above boolean vector into a 1440×720 Matrix{Bool}:

inW1 = permutedims(reshape(inW, length(lat), length(lon)))
See comparison herein with MWE :
using PolygonOps, StaticArrays

lonW = [50, 100, 100, 50., 50.]
latW = [0, 0, 30, 30., 0.]
polygon = SVector.(lonW, latW)    # boundary of the polygon

lon  = 20:0.25:379.75;   # 1440 values
lat  = -90:0.25:89.75;   #  720 values
points = vec(SVector.(lon',lat))

inW = [inpolygon(p, polygon; in=true, on=false, out=false) for p in points]
inW1 = permutedims(reshape(inW, length(lat), length(lon)))

inW2 = [inpolygon((x, y), polygon; in=true, on=false, out=false) for x in lon, y in lat]
inW1 == inW2  # true

using Plots; gr()
plot(heatmap(lon, lat,inW1'), heatmap(lon,lat, inW2'), cbar=false)
1 Like

ndgrid in DIVAnd is just provided as a utility tool: DIVAnd works on a possibly curvilinear grid (for which you need to have the information of the position of each grid point). To help users create a grid in case they want to use a uniform rectangular grid, we provide ndgrid to do that for them.

1 Like