Indexing a 4-D array using a logical matrix

Hello, I have a few 4-D arrays, originally netCDF files, which represent mapping data. Each 4-D array, or map, is lat x lon x percentile x frequency. I have the actual lat, lon, percentile, and frequency values sorted as seperate vectors.

Each map is actually 1144 x 1051 x 7 x 3, so for practice:

map=rand(1144,1051,7,3)

I have a polygon I want to subset the data by, and I have created an index by:

eez=shaperead('./ShapeFiles/eez/eez.shp');
lat = ncread(ncfile, '/lat');
lon = ncread(ncfile, '/lon');
[LonGrid, LatGrid] = meshgrid(lon, lat);  %make a grid with the latlongs
Y=eez.Y(~isnan(eez.Y))'; %removing NAs and transposing the array
X=eez.X(~isnan(eez.X))';
inPoly = inpolygon(LonGrid, LatGrid,X,Y);

How do I use the inPoly logical to subset the 4-D map, just for the lat/lon, but apply it equally to all other dimensions (percentile/frequency), so I have a resulting 4-D array with the subsetted data? If I just index it by:

test=map(inPoly);

It returns an array that I am unsure how to apply to the original 4-D matrix. Ideally, it would return a 4D matrix with NaNs outside of the polygon.

Thanks in advance!

[Matt J]

How do I use the inPoly logical to subset the 4-D map, just for the lat/lon, but apply it equally to all other dimensions (percentile/frequency), so I have a resulting 4-D array with the subsetted data

It is not clear why you would expect a 4-D array, rather than a 3-D array, as a result. When you extract the interior pixels of an arbitrary polygon, the pixels do not have any natural rectangular shape, and must be organized as a 1D list.

[Emily T. Griffiths]

Good point! I will edit the question. I would want there to be NaNs outside the polygon.

Answer 1

[Matt J]

Assuming you are trying to extract the polygon interiors, you could do,

[~,~,p,q]=size(map);

Map=reshape(map,[],p*q);

test=reshape( Map(inPoly,:) , p, q );

If you are only trying to mask the array so that everything outside the polygon is set to zero, then you could just do instead,

test=map.*inpPoly;

[Emily T. Griffiths]

Thanks! I think ideally I would like to set everything outside the polygon to NaN. However, 0 is an inprobably number in this dataset, so this may work, as I can then reassign everything that is 0 to NaN. :)

[Matt J]

Thanks!

You are quite welcome, but lease Accept-click the answer if your question is resolved.

ideally I would like to set everything outside the polygon to NaN

That can be done as below:

[~,~,p,q]=size(map);
Map=reshape(map,[],p*q);
Map(~inPoly,:)=nan;

map=reshape(Map, size(map));

or possibly also,

map=map./inPoly;
map(~isfinite(map))=nan;