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

**URL:** <https://discourse.julialang.org/t/how-to-check-if-my-gridded-x-and-y-are-inside-a-polygon/70728>\
**Category:** General Usage\
**Tags:** question, gmt, polygons\
**Created:** [November 1, 2021, 10:00am UTC](https://discourse.julialang.org/t/how-to-check-if-my-gridded-x-and-y-are-inside-a-polygon/70728 "2021-11-01T10:00:21Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![leon](https://avatars.discourse-cdn.com/v4/letter/l/dc4da7/32.png) [@leon](https://discourse.julialang.org/u/leon)\
**Post date:** [November 1, 2021, 10:00am UTC](https://discourse.julialang.org/t/how-to-check-if-my-gridded-x-and-y-are-inside-a-polygon/70728/1 "2021-11-01T10:00:21Z")

</div>

I have gridded X and Y that are defined by

```julia
lon = 20:0.25:379.75;
lat = -90:0.25:89.75;
X, Y = ndgrid(lon, lat).

```

My goal is to create a Boolean array (named `Index`) the same size as X (or Y) and their values should be 1 if they are inside a polygon. That way I can apply it to my Variable grid like: Variable[Index] .= NaN;

Below is what I found from the internet, but the generated Boolean array is a column data, instead of the size of X.

```julia
using PolygonOps
using StaticArrays

polygon = SVector.(lonW, latW); # boundary of the polygon
lon = 20:0.25:379.75;
lat = -90:0.25:89.75;
points = vec(SVector.(lon',lat));

inW = [inpolygon(p, polygon; in=true, on=false, out=false) for p in points];

```

How can I replace the `points` above with my gridded X and Y, so as to achieve my above goal? Thanks!

---

<div class="post-metadata">

**Author:** ![gustaphe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gustaphe/32/18174_2.png) [@gustaphe](https://discourse.julialang.org/u/gustaphe)\
**Post date:** [November 1, 2021, 10:26am UTC](https://discourse.julialang.org/t/how-to-check-if-my-gridded-x-and-y-are-inside-a-polygon/70728/2 "2021-11-01T10:26:29Z")

</div>

What is `ndgrid`?

---

<div class="post-metadata">

**Author:** ![gustaphe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gustaphe/32/18174_2.png) [@gustaphe](https://discourse.julialang.org/u/gustaphe)\
**Post date:** [November 1, 2021, 10:42am UTC](https://discourse.julialang.org/t/how-to-check-if-my-gridded-x-and-y-are-inside-a-polygon/70728/3 "2021-11-01T10:42:36Z")

</div>

Not knowing what `ndgrid` does (or more importantly, where it’s defined), I think

```julia
inW = [inpolygon((x, y), polygon; in=true, on=false, out=false) for x in lon, y in lat]

```

should do what you want. In general, I haven’t seen very many strong use cases for “grids” in Julia, broadcasting and loops are just that much better.

---

<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:** [November 1, 2021, 11:38am UTC](https://discourse.julialang.org/t/how-to-check-if-my-gridded-x-and-y-are-inside-a-polygon/70728/4 "2021-11-01T11:38:52Z")

</div>

Let’s dig a square hole on Atlantic

```julia
using GMT

# Load the earth grid at 15 arc minutes
G = gmtread("@earth_relief_15m");

# A polygon
pol = [-58 13; -58 37; -20 37; -20 13; -58 13];

# Compute the mask. Here the ``inc`` and ``registration`` should be 
# automatically assigned via the ``G`` argument in ``region`` but some bug ...
# Anyway, understanding *pixel vs grid* registration is fundamental
# https://docs.generic-mapping-tools.org/dev/cookbook/options.html#pixel-registration
mask = grdmask(pol, region=G, inc=G.inc, out_edge_in=(1,1,NaN), registration=:p);

# Now let's dig that hole
holed_earth = G * mask;

imshow(holed_earth, proj=:guess, shade=true)

```

 ![GMTjl_tmp-fs8](https://global.discourse-cdn.com/julialang/original/3X/f/0/f0256683a2f54af9164d767f5166ba92833ea29e.jpeg)

---

<div class="post-metadata">

**Author:** ![leon](https://avatars.discourse-cdn.com/v4/letter/l/dc4da7/32.png) [@leon](https://discourse.julialang.org/u/leon)\
**Post date:** [November 1, 2021, 11:55am UTC](https://discourse.julialang.org/t/how-to-check-if-my-gridded-x-and-y-are-inside-a-polygon/70728/5 "2021-11-01T11:55:36Z")

</div>

Many thanks, All.

`ndgrid` is similar to `meshgrid`. It is a function within the DIVAnd.jl package.

---

<div class="post-metadata">

**Author:** ![leon](https://avatars.discourse-cdn.com/v4/letter/l/dc4da7/32.png) [@leon](https://discourse.julialang.org/u/leon)\
**Post date:** [November 1, 2021, 12:00pm UTC](https://discourse.julialang.org/t/how-to-check-if-my-gridded-x-and-y-are-inside-a-polygon/70728/6 "2021-11-01T12:00:02Z")

</div>

Many thanks!

Is `G` a n1 x n2 gridded data? Can I apply the same thing to another grid that is produced with a different set of longitude and latitude? Or this is only applicable to internal GMT stored grids?

---

<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:** [November 1, 2021, 12:00pm UTC](https://discourse.julialang.org/t/how-to-check-if-my-gridded-x-and-y-are-inside-a-polygon/70728/7 "2021-11-01T12:00:07Z")

</div>

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

---

<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:** [November 1, 2021, 12:03pm UTC](https://discourse.julialang.org/t/how-to-check-if-my-gridded-x-and-y-are-inside-a-polygon/70728/8 "2021-11-01T12:03:51Z")

</div>

> [@leon](#):
>
> Or this is only applicable to internal GMT stored grids?

_GMT stored grids_ is just a Julia type. I already pointed you to [them](https://www.generic-mapping-tools.org/GMT.jl/dev/types/#Grid-type). It has a `G.z` member that is a plain 2D array.

---

<div class="post-metadata">

**Author:** ![leon](https://avatars.discourse-cdn.com/v4/letter/l/dc4da7/32.png) [@leon](https://discourse.julialang.org/u/leon)\
**Post date:** [November 1, 2021, 12:28pm UTC](https://discourse.julialang.org/t/how-to-check-if-my-gridded-x-and-y-are-inside-a-polygon/70728/9 "2021-11-01T12:28:38Z")

</div>

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.

---

<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:** [November 1, 2021, 12:40pm UTC](https://discourse.julialang.org/t/how-to-check-if-my-gridded-x-and-y-are-inside-a-polygon/70728/10 "2021-11-01T12:40:44Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![gustaphe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gustaphe/32/18174_2.png) [@gustaphe](https://discourse.julialang.org/u/gustaphe)\
**Post date:** [November 1, 2021, 1:10pm UTC](https://discourse.julialang.org/t/how-to-check-if-my-gridded-x-and-y-are-inside-a-polygon/70728/11 "2021-11-01T13:10:30Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [November 1, 2021, 2:35pm UTC](https://discourse.julialang.org/t/how-to-check-if-my-gridded-x-and-y-are-inside-a-polygon/70728/12 "2021-11-01T14:35:36Z")

</div>

> [@gustaphe](#):
>
> In general, I haven’t seen very many strong use cases for “grids” in Julia

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?

---

<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:** [November 1, 2021, 2:57pm UTC](https://discourse.julialang.org/t/how-to-check-if-my-gridded-x-and-y-are-inside-a-polygon/70728/13 "2021-11-01T14:57:14Z")

</div>

And I bet your [gorgeous golden stripes](https://discourse.julialang.org/t/3d-heatmap/70658/4), too.

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [November 1, 2021, 3:08pm UTC](https://discourse.julialang.org/t/how-to-check-if-my-gridded-x-and-y-are-inside-a-polygon/70728/14 "2021-11-01T15:08:43Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![leon](https://avatars.discourse-cdn.com/v4/letter/l/dc4da7/32.png) [@leon](https://discourse.julialang.org/u/leon)\
**Post date:** [November 1, 2021, 3:15pm UTC](https://discourse.julialang.org/t/how-to-check-if-my-gridded-x-and-y-are-inside-a-polygon/70728/15 "2021-11-01T15:15:09Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![gustaphe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gustaphe/32/18174_2.png) [@gustaphe](https://discourse.julialang.org/u/gustaphe)\
**Post date:** [November 1, 2021, 3:31pm UTC](https://discourse.julialang.org/t/how-to-check-if-my-gridded-x-and-y-are-inside-a-polygon/70728/16 "2021-11-01T15:31:22Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [November 1, 2021, 3:37pm UTC](https://discourse.julialang.org/t/how-to-check-if-my-gridded-x-and-y-are-inside-a-polygon/70728/17 "2021-11-01T15:37:05Z")

</div>

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

---

<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:** [November 1, 2021, 3:55pm UTC](https://discourse.julialang.org/t/how-to-check-if-my-gridded-x-and-y-are-inside-a-polygon/70728/18 "2021-11-01T15:55:18Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![aramirezreyes](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/aramirezreyes/32/42573_2.png) [@aramirezreyes](https://discourse.julialang.org/u/aramirezreyes)\
**Post date:** [November 1, 2021, 4:24pm UTC](https://discourse.julialang.org/t/how-to-check-if-my-gridded-x-and-y-are-inside-a-polygon/70728/19 "2021-11-01T16:24:37Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![gustaphe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gustaphe/32/18174_2.png) [@gustaphe](https://discourse.julialang.org/u/gustaphe)\
**Post date:** [November 1, 2021, 5:54pm UTC](https://discourse.julialang.org/t/how-to-check-if-my-gridded-x-and-y-are-inside-a-polygon/70728/20 "2021-11-01T17:54:41Z")

</div>

Not really.

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

[Next page](https://discourse.julialang.org/t/how-to-check-if-my-gridded-x-and-y-are-inside-a-polygon/70728.md?page=2)