BenWalleshauser
commented
2 years ago Hello,
Thank you for reaching out and voicing your concerns regarding the code! The current difficulty you are experiencing is due to how neighbors are found for indices at the top of the map.
Following the 120x240 matrix used in the paper, let's say we want to find the neighboring indices for the first pack:
We can first find the neighbors directly surrounding the pack (please note the indices in the East):
Now there is a challenge to find the neighbors for the indices on the top of the map, as on a globe they would all converge at the pole. But it would not be computationally efficient to include all of the points on the top of the map (total size of input for RC would be >240). Therefore the Neighbors.m file currently will select the corresponding index that is symmetric about the prime meridian (note the same computations were neglected in the Neighbors.m file for the South Pole due to the used dataset being entirely land/ice there).
Using a combination of the FindPackIndices.m and Neighbors.m file this is also what is produced (for a total of 34 indices).
clear
clc
nBox = 120; %Number of latitudes
mBox = 240; %Number of longitudes
%Number of rows and columns per RC (can change)
nPack = 4;
mPack = 4;
NumPacks = (nBox*mBox)/(nPack*mPack); %Number of packs, each pack will have its own RC
NaNset = 250;
%function PackInd = FindPackIndices(Avgs,nBox,mBox,nPack,mPack,NaNset)
%Returns the matrix indices that are in each pack. Each column of indices
%contains the points within a pack.
if nBox/nPack ~= ceil(nBox/nPack) | mBox/mPack ~= ceil(mBox/mPack)
error('Number of Boxes needs to be divisible by num of RCs (Packs)')
end
%Finding the indices of the boxes per pack/RC
PackInd = zeros(nPack*mPack,NumPacks);
d = 1;
c = 0;
b = 0;
e = 0;
mult = 0;
mult2 = 0;
for i = 1:NumPacks
for j = 1:mPack
for k = 1:nPack
%Cycling through each pack:
PackInd(b+1,i) = (j-1)*nBox + k + mult2*nPack + mult*nBox*mPack;
b = b+1;
%if SST_total(1,PackInd(b,i)) == NaNset %This just disregards points on land from the creation of the pack
%PackInd(b,i) = 0;
%end
end
if b >= nPack*mPack
mult2 = mult2+1;
b = 0;
end
end
c = c + 1;
if c >= nBox/nPack
mult = mult+1;
c = 0;
mult2 = 0;
end
end
%end
%% Finding Neighbors
NeighborInd = zeros((nPack+2)*(mPack+2),NumPacks);
%Finding of neighboring indices
for p = 1:NumPacks %creating RC for each group
iPackInd = nonzeros(PackInd(:,p)); %The indices that are going to be predicted with the given RC
CloudInd = zeros(8,length(iPackInd)); %The indices that surround the pack
%Finding the indices that are going to be considered in RC
for z = 1:length(iPackInd)
nbInds = Neighbors(iPackInd(z),nBox,mBox); %The neighboring indices for a given point
%for i = 1:length(nbInds) %Getting rid of neighbor if it is on land
%if SST_total(1,nbInds(i)) == NaNset
%nbInds(i) = 0;
%end
%end
CloudInd(1:length(nbInds),z) = nbInds; %All of the neighboring indices
end
CloudInd = reshape(CloudInd,[],1);
TotalInd = [iPackInd; CloudInd];
TotalInd = unique(TotalInd,'stable'); %Getting rid of repeat Neighbors while preserving order of indices
TotalInd = TotalInd(TotalInd~=0); %Total Indices being condsidered in the given RC (nonzero)
NeighborInd(1:length(TotalInd),p) = TotalInd;
end
PackOne = sort(nonzeros(NeighborInd(:,1))); %Total indices in first pack
%% Neighbors function as in current code
function [neighbors] = Neighbors(ind, nBox, mBox)
%Returns the neighboring indices surrounding a point on the globe. Disregards the
%neighbors that wrap around on the bottom of the map as they are all land
%values (due to Antarctica).
k = 0; %index of neighbors, increases as conditions met
%% Top, Bottom, Left, and Right
%Neighbor beneath
Bottom = 0;
if ind/nBox ~= ceil(ind/nBox) %ie not at bottom of map
k = k+1;
Bottom = 1;
neighbors(k) = ind+1;
end
%Neighbor above
Top = 0;
if (ind-1)/nBox ~= ceil((ind-1)/nBox) & ind ~= 1 %If point is not on top of the map
k = k+1;
Top = 1;
neighbors(k) = ind-1;
end
TopMap = 0;
if mod(ind-1,nBox) == 0 %If point is on top of the map
k = k+1;
TopMap = 1;
myCol = (ind-1)/nBox+1;
neighborCol = 240-myCol+1;
neighbors(k) = (neighborCol-1)*nBox+1;
end
%Neighbor right
if ind <= nBox*mBox-nBox %If point is not on the right edge of map
k = k+1;
neighbors(k) = ind + nBox;
end
if ind > nBox*mBox-nBox %If point is on the right edge of map
k = k+1;
neighbors(k) = ind - (nBox*mBox-nBox);
end
%Neighbor left
if ind >= nBox+1 %If point is not on the left edge of the map
k = k+1;
neighbors(k) = ind - nBox;
end
if ind < nBox+1 %If point is on the left edge of the map
k = k+1;
neighbors(k) = ind + (nBox*mBox-nBox);
end
%% Corners
%Top Left
if Top == 1 %If not on top edge
k = k+1;
neighbors(k) = ind-nBox-1;
if neighbors(k) < 1 %If on left edge therefore the neighbor wraps around
neighbors(k) = neighbors(k) + nBox*mBox;
end
end
if TopMap == 1; %If on top edge
neighborCol = 240-myCol+1+1;
if neighborCol ~= mBox+1 %Would be it's own neighbor if condition true
k = k+1;
neighbors(k) = (neighborCol-1)*nBox+1;
end
end
%Top Right
if Top == 1 %If not on top edge
k = k+1;
neighbors(k) = ind+nBox-1;
if neighbors(k) > nBox*mBox %If on right edge therefore the neighbor wraps around
neighbors(k) = neighbors(k) - nBox*mBox;
end
end
if TopMap == 1; %If on top edge
neighborCol = 240-myCol+1-1;
if neighborCol ~= 0 %Would be it's own neighbor if condition true
k = k+1;
neighbors(k) = (neighborCol-1)*nBox+1;
end
end
%Bottom Right
if Bottom == 1
k = k+1;
neighbors(k) = ind+nBox+1;
if neighbors(k) > nBox*mBox %If on right edge therefore the neighbor wraps around
neighbors(k) = neighbors(k) - nBox*mBox;
end
end
%Bottom Left
if Bottom == 1
k = k+1;
neighbors(k) = ind-nBox+1;
if neighbors(k) < 1 %If on left edge therefore the neighbor wraps around
neighbors(k) = neighbors(k) + nBox*mBox;
end
end
end
If you'd like to disregard the neighbors found across the pole, please use this revised version of Neighbors.m:
function [neighbors] = Neighbors(ind, nBox, mBox)
%Returns the neighboring indices surrounding a point on the globe. Disregards the
%neighbors that wrap around on the top and bottom of the map.
k = 0; %index of neighbors, increases as conditions met
%% Top, Bottom, Left, and Right
%Neighbor beneath
Bottom = 0;
if ind/nBox ~= ceil(ind/nBox) %ie not at bottom of map
k = k+1;
Bottom = 1;
neighbors(k) = ind+1;
end
%Neighbor above
Top = 0;
if (ind-1)/nBox ~= ceil((ind-1)/nBox) & ind ~= 1 %If point is not on top of the map
k = k+1;
Top = 1;
neighbors(k) = ind-1;
end
%Neighbor right
if ind <= nBox*mBox-nBox %If point is not on the right edge of map
k = k+1;
neighbors(k) = ind + nBox;
end
if ind > nBox*mBox-nBox %If point is on the right edge of map
k = k+1;
neighbors(k) = ind - (nBox*mBox-nBox);
end
%Neighbor left
if ind >= nBox+1 %If point is not on the left edge of the map
k = k+1;
neighbors(k) = ind - nBox;
end
if ind < nBox+1 %If point is on the left edge of the map
k = k+1;
neighbors(k) = ind + (nBox*mBox-nBox);
end
%% Corners
%Top Left
if Top == 1 %If not on top edge
k = k+1;
neighbors(k) = ind-nBox-1;
if neighbors(k) < 1 %If on left edge therefore the neighbor wraps around
neighbors(k) = neighbors(k) + nBox*mBox;
end
end
%Top Right
if Top == 1 %If not on top edge
k = k+1;
neighbors(k) = ind+nBox-1;
if neighbors(k) > nBox*mBox %If on right edge therefore the neighbor wraps around
neighbors(k) = neighbors(k) - nBox*mBox;
end
end
%Bottom Right
if Bottom == 1
k = k+1;
neighbors(k) = ind+nBox+1;
if neighbors(k) > nBox*mBox %If on right edge therefore the neighbor wraps around
neighbors(k) = neighbors(k) - nBox*mBox;
end
end
%Bottom Left
if Bottom == 1
k = k+1;
neighbors(k) = ind-nBox+1;
if neighbors(k) < 1 %If on left edge therefore the neighbor wraps around
neighbors(k) = neighbors(k) + nBox*mBox;
end
end
endPlease let me know if you still experience issues!
Best regards, Ben W.
renierts
Hello,
thanks for sharing the code for your paper. I have a question about the computation of the neighborhood of a patch. To debug this, I created the following matrix containing all of the indices (regardless the current point is land or water):
I can pad this matrix in order to contain surrounding neighbors:
Now, I would like to follow you and extract 4x4 patches together with the neigbors. This breaks down to extract a 6x6 submatrix, shift by 4 and extract the next 6x6 matrix. The first patch with its neighbors would look like the following:
Two issues are:
Could you please help me to identify potential mistakes in my way of thinking?