Home > atmlab > geoplots > gmt > polygoninize_regions.m

polygoninize_regions

PURPOSE ^

POLYGONINIZE_REGIONS Turn adjacent regions (boxes) into polygon for plotting

SYNOPSIS ^

function pspoly = polygoninize_regions(boxes)

DESCRIPTION ^

 POLYGONINIZE_REGIONS Turn adjacent regions (boxes) into polygon for plotting

 PURPOSE: To set up for plotting many regions defined as box coordinates. This
          function make sure that extra lines aren't drawn right next to each
          other for essentially the same region. This is usefull if you have a
          complicated region made of of several boxes.
          e.g. in eliasson11_assessing_acp
         

 IN      boxes  [%f %s]: (array containing coordinates for one or several regions)
         e.g. boxes = [box1; box2; box3], where box1 has the form 
         [bottom left corner, top right corner];

 OUT     struct      Structure of arguments to be used for plotting the
                     regions  on a map

 Created by Salomon Eliasson
 $Id: polygoninize_regions.m 6862 2011-04-17 20:27:55Z seliasson $

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

DOWNLOAD ^

polygoninize_regions.m

SOURCE CODE ^

0001 function pspoly = polygoninize_regions(boxes)
0002 % POLYGONINIZE_REGIONS Turn adjacent regions (boxes) into polygon for plotting
0003 %
0004 % PURPOSE: To set up for plotting many regions defined as box coordinates. This
0005 %          function make sure that extra lines aren't drawn right next to each
0006 %          other for essentially the same region. This is usefull if you have a
0007 %          complicated region made of of several boxes.
0008 %          e.g. in eliasson11_assessing_acp
0009 %
0010 %
0011 % IN      boxes  [%f %s]: (array containing coordinates for one or several regions)
0012 %         e.g. boxes = [box1; box2; box3], where box1 has the form
0013 %         [bottom left corner, top right corner];
0014 %
0015 % OUT     struct      Structure of arguments to be used for plotting the
0016 %                     regions  on a map
0017 %
0018 % Created by Salomon Eliasson
0019 % $Id: polygoninize_regions.m 6862 2011-04-17 20:27:55Z seliasson $
0020 
0021 %% Make a grid of logicals
0022 res = .1;
0023 tmplon = -180-res:res:180+res;
0024 tmplat = -90-res:res:90+res;
0025 tmpgrid = false(length(tmplat),length(tmplon));
0026 [lon,lat]=meshgrid(tmplon(2:end-1),tmplat((2:end-1)));
0027 
0028 % cond1 within the regions
0029 for i=1:size(boxes,1)
0030     latindex= tmplat >= boxes(i,1) & tmplat <= boxes(i,3);
0031     lonindex= tmplon >= boxes(i,2) & tmplon <= boxes(i,4);
0032     tmpgrid(latindex,lonindex)=true;
0033 end
0034 
0035 j = 2:length(tmplat)-1;
0036 i = 2:length(tmplon)-1;
0037 
0038 %% grid
0039 %|------|------|------| |------|------|------|
0040 %| j-1, | j-1, | j-1, | |  1   | 2    | 3    |
0041 %|   i-1|   i  |   i+1| |      |      |      |
0042 %|------|------|------| |------|------|------|
0043 %| j,   | j,i  | j,   | |  4   | i,j  |  5   |
0044 %|  i-1 |      |  i+1 | |      |      |      |
0045 %|------|------|------| |------|------|------|
0046 %| j+1, | j+1, | j+1, | |  6   |  7   |  8   |
0047 %|   i-1|   i  |   i+1| |      |      |      |
0048 %|------|------|------| |------|------|------|
0049 
0050 cond1= tmpgrid(j,i) & tmpgrid(j-1,i-1);
0051 cond2= tmpgrid(j,i) & tmpgrid(j-1,i);
0052 cond3= tmpgrid(j,i) & tmpgrid(j-1,i+1);
0053 cond4= tmpgrid(j,i) & tmpgrid(j,i-1);
0054 cond5= tmpgrid(j,i) & tmpgrid(j,i+1);
0055 cond6= tmpgrid(j,i) & tmpgrid(j+1,i-1);
0056 cond7= tmpgrid(j,i) & tmpgrid(j+1,i);
0057 cond8= tmpgrid(j,i) & tmpgrid(j+1,i+1);
0058 
0059 %% CORNERS
0060 % remember grid is upside down
0061 tl = tmpgrid(2:end-1,2:end-1) & ~(cond4 | cond2); %bl
0062 tr = tmpgrid(2:end-1,2:end-1) & ~(cond5 | cond2); %br
0063 bl = tmpgrid(2:end-1,2:end-1) & ~(cond4 | cond7); %tl
0064 br = tmpgrid(2:end-1,2:end-1) & ~(cond5 | cond7); %tr
0065 itl = xor(cond1,cond4 & cond2) & ~tmpgrid(j-1,i-1); % ibl
0066 itr = xor(cond3,cond5 & cond2) & ~tmpgrid(j-1,i+1); % ibr
0067 ibl = xor(cond6,cond4 & cond7) & ~tmpgrid(j+1,i-1); % itl
0068 ibr = xor(cond8,cond5 & cond7) & ~tmpgrid(j+1,i+1); % itr
0069 
0070 %% TRACING
0071 corners = tl | tr | bl | br | itl | itr | ibl | ibr;
0072 [rows,cols]=find(corners);
0073 points = [lat(corners),lon(corners),rows,cols];
0074 
0075 tmpgrid = tmpgrid(2:end-1,2:end-1);
0076 
0077 di=-1;
0078 j = 1;
0079 curpoints = [points,false(size(cols))];
0080 psxy = cell(10,1);
0081 while 1
0082     if isempty(curpoints)
0083         break
0084     end
0085     i = 1;
0086     psxy{j} = [curpoints(i,2),curpoints(i,1)];
0087     curpoints(i,5)=true;
0088     while  1
0089         di = finddirection(tmpgrid,curpoints,i,di);
0090         i = nextcoord(di,curpoints,i);
0091         if curpoints(i,5)
0092             break
0093         end
0094         psxy{j} = [psxy{j}; [curpoints(i,2),curpoints(i,1)]];
0095         curpoints(i,5)=true;
0096     end
0097     curpoints = curpoints(~curpoints(:,5),:);
0098     j = j+1;
0099 end
0100 
0101 pspoly = psxy(1:j-1); %pspoly is the name used in gmt_plot
0102 
0103 %% SUBFUNCTIONS
0104 %     |||||||
0105 %     vvvvvvv
0106 
0107 function d = finddirection(tmpgrid,points,i,dir)
0108 %% finddirection
0109 
0110 [d,f] = specialcase(tmpgrid,points,i,dir);
0111 
0112 if f, return, end
0113 
0114 for j=1:4
0115     d = mod(d+1,4);
0116     if d == dir
0117         break
0118     end
0119     if dir~=-1 && mod(dir+2,4)==d, continue, end
0120     switch d
0121         case 0 % west, left in matrix
0122             beside=~isempty(points(points(:,3) == points(i,3) & points(:,4) == points(i,3)-1,:));
0123             west = sum(tmpgrid(points(i,3)-1:points(i,3)+1,points(i,4)-1))==2;
0124             if beside || west
0125                 break
0126             end
0127         case 1 %north, down in matrix
0128             beside=~isempty(points(points(:,4) == points(i,4) & points(:,3) == points(i,3)+1,:));
0129             north = sum(tmpgrid(points(i,3)+1,points(i,4)-1:points(i,4)+1))==2;
0130             if beside || north
0131                 break
0132             end
0133         case 2 %east, right in matrix
0134             beside=~isempty(points(points(:,3) == points(i,3) & points(:,4) == points(i,3)+1,:));
0135             east = sum(tmpgrid(points(i,3)-1:points(i,3)+1,points(i,4)+1))==2;
0136             if beside || east
0137                 break
0138             end
0139         case 3 %south, up in matrix
0140             beside=~isempty(points(points(:,4) == points(i,4) & points(:,3) == points(i,3)-1,:));
0141             south = sum(tmpgrid(points(i,3)-1,points(i,4)-1:points(i,4)+1))==2;
0142             if beside || south
0143                 break
0144             end
0145     end
0146 end
0147 
0148 
0149 function a = nextcoord(di,points,i)
0150 %% nextcoord
0151 
0152 ref = points(i,3:4);
0153 
0154 k = 1:length(points);
0155 switch di
0156     case 0 % west, left in matrix
0157         index = ref(1,1)==points(:,3)&ref(1,2)>points(:,4);
0158     case 1 %north, down in matrix
0159         index = ref(1,2)==points(:,4)&ref(1,1)<points(:,3);
0160     case 2 %east, right in matrix
0161         index = ref(1,1)==points(:,3)&ref(1,2)<points(:,4);
0162     case 3 %south, up in matrix
0163         index = ref(1,2)==points(:,4)&ref(1,1)>points(:,3);
0164 end
0165 
0166 k = k(index);
0167 test = sum(repmat(ref,sum(index),1)-points(index,3:4),2);
0168 a = k(abs(test)==min(abs(test)));
0169 if isempty(a)
0170     error('gmtlab:error','isempty(a)')
0171 end
0172 
0173 function [dir,f] = specialcase(tmpgrid,points,i,dir)                            
0174 %% SPECIAL CASE
0175 
0176 a = points(i,4) == 1;              % you are on the western edge of the map
0177 b = points(i,4) == size(tmpgrid,2);% you are on the eastern edge of the map
0178 c = points(i,3) == 1;              % you are on the southern edge of the map
0179 d = points(i,3) == size(tmpgrid,1);% you are on the northern edge of the map
0180 f = false;
0181 if a || b || c || d
0182     if a
0183         % north or east
0184         if checknorth(tmpgrid,points,i)
0185             dir = 1;
0186         else dir = 2;
0187         end
0188     end
0189     if b
0190         % south or west
0191         if checksouth(tmpgrid,points,i)
0192             dir = 3;
0193         else dir = 0;
0194         end
0195     end
0196     if c
0197         % west or north
0198         if checkwest(tmpgrid,points,i)
0199             dir = 0;
0200         else dir = 1;
0201         end
0202     end
0203     if d
0204         % east or south
0205         if checkeast(tmpgrid,points,i)
0206             dir = 2;
0207         else dir = 3;
0208         end
0209     end
0210     f = true;
0211 end
0212 
0213 
0214 function out = checknorth(tmpgrid,points,i)                                     
0215 %% CHECK NORTH
0216 
0217 out = false;
0218 ref = points(i,3:4);
0219 k = 1:length(points);
0220 
0221 % check if any points exist in this direction
0222 cond = points(:,4)==ref(1,2)&points(:,3)>ref(1,1)&~points(:,5);
0223 if sum(cond)==0, return, end
0224 
0225 subset = points(cond,3:4);
0226 k = k(cond);
0227 test = sum(repmat(ref,size(subset,1),1)-subset,2);
0228 a = k(abs(test)==min(abs(test)));
0229 
0230 %test if it's there (western edge)
0231 vect = tmpgrid(ref(1,1):points(a,3),1);
0232 out = sum(vect,1)==length(vect);
0233 
0234 function out = checksouth(tmpgrid,points,i)                                     
0235 %% CHECK SOUTH
0236 
0237 out=false;
0238 ref = points(i,3:4);
0239 k = 1:length(points);
0240 
0241 % check if any points exist in this direction
0242 cond = points(:,4)==ref(1,2)&points(:,3)<ref(1,1)&~points(:,5);
0243 if sum(cond)==0, return, end
0244 
0245 subset = points(cond,3:4);
0246 k = k(cond);
0247 test = sum(repmat(ref,size(subset,1),1)-subset,2);
0248 a = k(abs(test)==min(abs(test)));
0249 
0250 %test if it's there (eastern edge)
0251 vect = tmpgrid(ref(1,1):points(a,3),end);
0252 out = sum(vect,1)==length(vect);
0253 
0254 function out = checkwest(tmpgrid,points,i)                                      
0255 %% CHECK WEST
0256 
0257 out=false;
0258 ref = points(i,3:4);
0259 k = 1:length(points);
0260 
0261 cond = points(:,3)==ref(1,1)&points(:,4)<ref(1,2)&~points(:,5);
0262 if sum(cond)==0, return, end
0263 
0264 subset = points(cond,3:4);
0265 k = k(cond);
0266 test = sum(repmat(ref,size(subset,1),1)-subset,2);
0267 a = k(abs(test)==min(abs(test)));
0268 
0269 %test if it's there (southern edge)
0270 vect = tmpgrid(1,ref(1,2):points(a,4));
0271 out = sum(vect,1)==length(vect);
0272 
0273 function out = checkeast(tmpgrid,points,i)                                      
0274 %% CHECK EAST
0275 
0276 out=false;
0277 ref = points(i,3:4);
0278 k = 1:length(points);
0279 
0280 cond = points(:,3)==ref(1,1)&points(:,4)>ref(1,2)&~points(:,5);
0281 if sum(cond)==0, return, end
0282 
0283 subset = points(cond,3:4);
0284 k = k(cond);
0285 test = sum(repmat(ref,size(subset,1),1)-subset,2);
0286 a = k(abs(test)==min(abs(test)));
0287 
0288 %test if it's there (northern edge)
0289 vect = tmpgrid(end,ref(1,2):points(a,4));
0290 out = sum(vect,1)==length(vect);

Generated on Mon 15-Sep-2014 13:31:28 by m2html © 2005