# 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:** 6\
**Page:** 2

<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, 6:07pm UTC](https://discourse.julialang.org/t/how-to-check-if-my-gridded-x-and-y-are-inside-a-polygon/70728/21 "2021-11-01T18:07:56Z")

</div>

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.

---

<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, 6:10pm UTC](https://discourse.julialang.org/t/how-to-check-if-my-gridded-x-and-y-are-inside-a-polygon/70728/22 "2021-11-01T18:10:32Z")

</div>

> [@gustaphe](#):
>
> Of course you need all of your samples, but you really don’t need to gridify `x` and `y` .

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)

---

<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, 7:15pm UTC](https://discourse.julialang.org/t/how-to-check-if-my-gridded-x-and-y-are-inside-a-polygon/70728/23 "2021-11-01T19:15:42Z")

</div>

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

---

<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, 7:31pm UTC](https://discourse.julialang.org/t/how-to-check-if-my-gridded-x-and-y-are-inside-a-polygon/70728/24 "2021-11-01T19:31:29Z")

</div>

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! 👍

---

<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, 8:55pm UTC](https://discourse.julialang.org/t/how-to-check-if-my-gridded-x-and-y-are-inside-a-polygon/70728/25 "2021-11-01T20:55:05Z")

</div>

> [@leon](#):
>
> `inW = [inpolygon(p, polygon; in=true, on=false, out=false) for p in points]`

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

```julia
inW1 = permutedims(reshape(inW, length(lat), length(lon)))

```

> **See comparison herein with MWE :**
>
> ```julia
> 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)
> 
> ```

---

<div class="post-metadata">

**Author:** ![JM\_Beckers](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jm_beckers/32/22482_2.png) [@JM\_Beckers](https://discourse.julialang.org/u/JM_Beckers)\
**Post date:** [November 3, 2021, 12:52pm UTC](https://discourse.julialang.org/t/how-to-check-if-my-gridded-x-and-y-are-inside-a-polygon/70728/26 "2021-11-03T12:52:08Z")

</div>

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.

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