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.
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!
Not knowing what ndgrid does (or more importantly, where it’s defined), I think
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.
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)
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?
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.
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.
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.
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.
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).
