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

General Usage
, ,
#7

you want to use GMT.jl and for other reasons as well) you should surface (or other GMT gridders) instead of ndgrid

1 Like
#8

GMT stored grids is just a Julia type. I already pointed you to them. It has a G.z member that is a plain 2D array.

1 Like
#9

Ccing @Alexander-Barth here.

I think ndgrid is the preferred gridding method by DIVAnd, am I right? I hope it will not be too much different from the surface function within GMT.

#10

I’m not comparing numeric merits (don’t even know that ndgrid, I thought it was the Matlab one) but you are having difficulties with the type conversions that would not exist if you computed the grid with a GMT method in first place.

#11

Next question: What is meshgrid? That’s not a Base function either.

If X, Y = ndgrid(x, y) are matrices such that X[i, j] == x[i] and Y[i, j] == y[j] (or vice versa), this is a very wasteful way to store data. For languages like Matlab it might make sense to create large matrices to represent repeated data like this, because they perform all their calculations matrixwise anyway. But Julia is intentionally designed to quickly and intelligently broadcast across the elements of vectors (or tuples, generators …). For instance, your lon and lat essentially store three numbers each: start, end and step. Turning them into a vector (lon_vec = collect(lon)) increases that number to length(lon) == 1438 and length(lat) == 720 respectively. That’s pretty wasteful. Now doing ndgrid spreads that same six-number (384 b) information over 2*length(lon)*length(lat) == 2070720 numbers (132 Mb). Unless of course the return type of ndgrid is just a wrapper of the original range, in which case I don’t see what use it is to begin with.

2 Likes
#12

I may be wrong, but from memory, some parametric 3D surface plots (using plotlyjs, pyplot, Makie, etc.) seem to require the creation of grids first?

#13

And I bet your gorgeous golden stripes, too.

1 Like
#14

@joa-quim, there may be job opportunities for a goldsmith and jewelry designer using Julia… :ring:

#15

I think your understanding is correct.

Basically ndgrid creates a matrix for both longitude and latitude. That way, the program knows the exact longitude and latitude values at each of the grid points.

That said, I agree the Julia way is more efficient.

#16

Plots.jl (which I’ve spent most of my time using) has plot(-2:0.1:2, 0:0.1:4, (x, y) -> exp(- x^2/4 - y^2/2)). Internally it might be using grids (depending on backend), which I would guess is at least not a bottleneck in the plotting flow, but at least you’re not required to drag the grids around as a user.

#17

Sure, but that is far from a parametric 3D surface like a torus or something fancy (and on top of which one may want to display some scalar field too).

#18

No need to go for fancy surfaces. If one has a measured quantity (e.g., temperature) on a regular mesh we must deal with grids.

1 Like
#19

And let’s not start talking about data on unstructured meshes!

1 Like
#20

Not really.

T = 273u"K" + 15u"K"*randn(501,501)
x = range(-5,5; length=501)u"cm"
y = range(-5,5; length=501)u"cm"

plot(x, y, T; st=:surface)

Of course you need all of your samples, but you really don’t need to gridify x and y.

#21

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.

#22

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)

#23

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

#24

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:

#25

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
#26

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