Home > atmlab > geographical > regionize.m

regionize

PURPOSE ^

% regionize

SYNOPSIS ^

function [outdata,outlat,outlon,weights] = regionize(data,lat,lon,corner,trim)

DESCRIPTION ^

% regionize

 PURPOSE: function that nan's out all data outside a given region bounds (or several regions)

 IN: data %2 or more dims starting with (lat,lon,...) or (lon,lat,...)
     lat  %vector
     lon  %vector
     corner = [blcorner,trcorner;blcorner,trcorner] % each corner is [lat,lon]
     trim   = true or false 
          (1=trim data and lat lons outside the region (default), 0=don't)

     NOTE REGARDING corner: If you want to use a predefined region, use
            e.g. corner = getPredefinedRegion('tropics') %look into function region list

     NOTE REGARDING trim: If you need the size of data matrix to remain unchanged, input trim=0

 OUT:    outdata
         outlat
         outlon
         weights % takes care of grids that are overlapping the edges
         flags = 

 USAGE: [data,lat,lon,weights,flags] = regionize(data,lat,lon,corner)
 
 NOTE: 1) If the data is gridded the spacing between lats and lon must be
          equidegree.
       2) Gridded data that is partly in and partly out of a region is
          included in outdata. Use the 3rd output argument to weight the
          data according to how much of the gridbox is in the region


 Created by Salomon Eliasson
 $Id: regionize.m 8740 2013-11-06 13:17:18Z seliasson $

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

DOWNLOAD ^

regionize.m

SOURCE CODE ^

0001 function [outdata,outlat,outlon,weights] = regionize(data,lat,lon,corner,trim)
0002 %% regionize
0003 %
0004 % PURPOSE: function that nan's out all data outside a given region bounds (or several regions)
0005 %
0006 % IN: data %2 or more dims starting with (lat,lon,...) or (lon,lat,...)
0007 %     lat  %vector
0008 %     lon  %vector
0009 %     corner = [blcorner,trcorner;blcorner,trcorner] % each corner is [lat,lon]
0010 %     trim   = true or false
0011 %          (1=trim data and lat lons outside the region (default), 0=don't)
0012 %
0013 %     NOTE REGARDING corner: If you want to use a predefined region, use
0014 %            e.g. corner = getPredefinedRegion('tropics') %look into function region list
0015 %
0016 %     NOTE REGARDING trim: If you need the size of data matrix to remain unchanged, input trim=0
0017 %
0018 % OUT:    outdata
0019 %         outlat
0020 %         outlon
0021 %         weights % takes care of grids that are overlapping the edges
0022 %         flags =
0023 %
0024 % USAGE: [data,lat,lon,weights,flags] = regionize(data,lat,lon,corner)
0025 %
0026 % NOTE: 1) If the data is gridded the spacing between lats and lon must be
0027 %          equidegree.
0028 %       2) Gridded data that is partly in and partly out of a region is
0029 %          included in outdata. Use the 3rd output argument to weight the
0030 %          data according to how much of the gridbox is in the region
0031 %
0032 %
0033 % Created by Salomon Eliasson
0034 % $Id: regionize.m 8740 2013-11-06 13:17:18Z seliasson $
0035 
0036 errID = 'atmlab:regionize:badInput';
0037 assert(nargin>=4,errID,'Not enought input arguments')
0038 
0039 % Default action is to trim away NaNs out side region
0040 if nargin==4, trim = true; end 
0041 
0042 [outdata,outlat,outlon,isgridded,flags] = setup(data,lat,lon);
0043 
0044 if isgridded
0045     % for the sake of grids overlapping the edges
0046     dlt = mean(diff(unique(lat)));
0047     dln = mean(diff(unique(lon)));
0048     
0049     
0050     sz = size(outdata);
0051 
0052     z.data = outdata(:,:,:);  %collapse to 3Doutdata;
0053 
0054     % NaN away everything that is not a region
0055     z.lat = repmat(outlat,1,length(outlon));
0056     z.lon = repmat(outlon,1,length(outlat))';
0057     
0058     [logical_field,periphery] = get_logicalfield(z,corner,dlt,dln);
0059     
0060     % Find weights for data on periphery of regions
0061     weights = findWeightsOnPeriphery(z.lat,z.lon,periphery,corner,dlt,dln);
0062     weights(logical_field) = 1;
0063 
0064     % Nan away everything outside region
0065     tmp =  permute(z.data,[3,1,2]); tsz = size(tmp);
0066     tmp = tmp(:,:);
0067     tmp(:,~logical_field&~periphery) = NaN;
0068     z.data = permute(reshape(tmp,tsz),[2,3,1]);
0069     
0070     %Trim away NaNs outside the region
0071     if trim
0072         [z.data,outlon,outlat,weights] = trimAwayNaNs_gridded(z.data,outlon,outlat,weights,corner,[dlt,dln]);
0073     end
0074     
0075     % Reshape the data array back to the original lat/lon format (first 2 dims)
0076     [outdata,outlat,outlon,weights] = useFlags(z.data,outlon,outlat,weights,flags);
0077     
0078     %resurrect matrix size to the original format
0079     outdata = reshape(outdata,[size(outdata,1),size(outdata,2),sz(3:end)]);
0080     
0081 else
0082     z.data        = outdata(:);
0083     z.lat         = outlat(:);
0084     z.lon         = outlon(:);
0085     logical_field = get_logicalfield(z,corner,0,0); %dlt and dln = 0
0086     
0087     outlat        = z.lat(logical_field);
0088     outlon        = z.lon(logical_field);
0089     outdata       = z.data(logical_field);
0090     
0091     if trim
0092        [outdata,outlon,outlat] = trimAwayNaNs_ungridded(outdata,outlon,outlat,corner); 
0093     end
0094 end
0095 
0096 %%%%%%%%%%%%%%%%%
0097 %% SUBFUNCTIONS
0098 %      |||||
0099 %      VVVVV
0100 function [outdata,outlat,outlon,ig,flags] = setup(data,lat,lon)
0101 %% setup
0102 
0103 errID = 'atmlab:reginize:badInput';
0104 
0105 a = size(data);
0106 ig = ~isequal(a,size(lat),size(lon));              % is gridded or not
0107 
0108 if ~ig
0109     outlon = lon(:);
0110     outlat = lat(:);
0111    outdata = data(:);
0112    flags = [];
0113    return
0114 end
0115 
0116 
0117 assert(isequal(a(1:2),[length(lat),length(lon)]) || isequal(a(1:2),[length(lon),length(lat)]),...
0118     errID,'Dimensions of data must be of form (lat,lon,...) or (lon,lat,...)')
0119 
0120 
0121 % make sure conventions are followed
0122 [flags,outlat,outlon,outdata] = standardize_geodata(lat(:),lon(:),data);
0123 
0124 function [logical_field,periphery] = get_logicalfield(in,corner,dlt,dln)
0125 %% get_logicalfield
0126 
0127 % fit exactly
0128 logical_field=false([size(in.data,1),size(in.data,2)]);
0129 for i = 1:size(corner,1)
0130     lt1 = in.lat >= corner(i,1) + dlt/2;
0131     lt2 = in.lat <= corner(i,3) - dlt/2;
0132     ln1 = in.lon >= corner(i,2) + dln/2;
0133     ln2 = in.lon <= corner(i,4) - dln/2;
0134 
0135     logical_field = logical_field | ( ln1&ln2&lt1&lt2 );
0136 end
0137 
0138 % on the periphery of the region
0139 periphery=false([size(in.data,1),size(in.data,2)]);
0140 for i = 1:size(corner,1)
0141     lt1 = in.lat >= corner(i,1) - dlt/2;
0142     lt2 = in.lat <= corner(i,3) + dlt/2;
0143     ln1 = in.lon >= corner(i,2) - dln/2;
0144     ln2 = in.lon <= corner(i,4) + dln/2;
0145 
0146     periphery = periphery | ( ln1&ln2&lt1&lt2 );
0147 end
0148 periphery(logical_field)=false;
0149 
0150 function weights = findWeightsOnPeriphery(lat,lon,periphery,C,dlt,dln)
0151 %% findWeightsOnPeriphery
0152 %
0153 % First NaNs away grids that are not on the perifery
0154 % Then loops over these points to find out how how they overlap the
0155 % region boundaries
0156 %
0157 % C = corners ([blcorner,trcorner;blcorner,trcorner]), dlt & dln are the gridsizes
0158 
0159 sz      = size(lat);
0160 weights = zeros(sz);
0161 num     = 1:length(lat(:));
0162 num     = reshape(num,sz);
0163 num     = num(periphery);
0164 
0165 % Now loop over remaining points. Loop also over the regions
0166 AREA = 0;
0167 for i = num'
0168     for j = 1:size(C,1)
0169         A = [lon(i)-dln/2,lat(i)-dlt/2,dln,dlt];
0170         B = [C(j,2),C(j,1),C(j,4)-C(j,2),C(j,3)-C(j,1)];
0171         AREA = AREA + rectint(A,B);
0172     end
0173     weights(i) = AREA/(dlt*dln); % normalized area. Max = 1;
0174     AREA = 0;
0175 end
0176 
0177 function [data,lon,lat,weights] = trimAwayNaNs_gridded(data,lon,lat,weights,corner,gsize)
0178 %% Trim dataset to get rid of NaNs outside the regions
0179 % If you don't need weights, make the 4th argument empty []
0180 % gsize is the size of the grid [dlat,dlon]
0181 
0182 gh = gsize/2; %half boxwidth
0183 
0184 % Use the corner input to trim away the data
0185 if ndims(data)==2
0186     data = data(lat>=min(corner(:,1))-gh(1) & lat<=max(corner(:,3))+gh(1),...
0187         lon>=min(corner(:,2))-gh(2) & lon<=max(corner(:,4))+gh(2));
0188 elseif ndims(data==3)
0189         data = data(lat>=min(corner(:,1))-gh(1) & lat<=max(corner(:,3))+gh(1),...
0190         lon>=min(corner(:,2))-gh(2) & lon<=max(corner(:,4))+gh(2),:);
0191 end
0192 
0193 if ~isempty(weights)
0194     weights = weights(lat>=min(corner(:,1))-gh(1) & lat<=max(corner(:,3))+gh(1),...
0195         lon>=min(corner(:,2)) -gh(2) & lon<=max(corner(:,4))+gh(2));
0196 end
0197 lat = lat(lat>=min(corner(:,1))-gh(1) & lat<=max(corner(:,3))+gh(1));
0198 lon = lon(lon>=min(corner(:,2))-gh(2) & lon<=max(corner(:,4))+gh(2));
0199 
0200 function [data,lon,lat,weights] = trimAwayNaNs_ungridded(data,lon,lat,corner)
0201 %% Trim dataset to get rid of NaNs outside the regions
0202 % If you don't need weights, make the 4th argument empty []
0203 % gsize is the size of the grid [dlat,dlon]
0204 
0205 % Use the corner input to trim away the data
0206 if ndims(data)==2
0207     data = data(lat>=min(corner(:,1)) & lat<=max(corner(:,3)) &...
0208         lon>=min(corner(:,2)) & lon<=max(corner(:,4)));
0209 elseif ndims(data==3)
0210         data = data(lat>=min(corner(:,1)) & lat<=max(corner(:,3)) &...
0211         lon>=min(corner(:,2)) & lon<=max(corner(:,4)),:);
0212 end
0213 
0214 lat = lat(lat>=min(corner(:,1)) & lat<=max(corner(:,3)));
0215 lon = lon(lon>=min(corner(:,2)) & lon<=max(corner(:,4)));
0216 
0217 function [data,lat,lon,weights] = useFlags(data,lon,lat,weights,flags)
0218 %% useFlags
0219 
0220 % if the orig lons are in the 0:360 regime. Put them BACK
0221 if flags.lon360
0222     lon = lon +(lon < 0)*360;
0223     [lon,lnindex] = sort(lon);
0224     data    = data(:,lnindex,:);
0225     weights = weights(:,lnindex);
0226 end
0227 %
0228 if flags.lon_descend
0229     [lon,lnindex]  = sort(lon,'descend');
0230     data = data(:,lnindex,:);
0231     weights = weights(:,lnindex);
0232 end
0233 
0234 % If the orig lats are DESCENDING. Put them BACK
0235 if flags.lat_descend
0236     [lat,ltindex]  = sort(lat,'descend');
0237     data = data(ltindex,:,:);
0238     weights = weights(ltindex,:);
0239 end
0240 
0241 % If orig data was data(lat,lon,...), Permute BACK
0242 if ~flags.permute
0243     data = permute(data,[2,1,3]);
0244     weights = permute(weights,[2,1]);
0245 end

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