All "regions" moved to REGIONS

csdms-contrib · Jul 1, 2016 · 41a2f13 · 41a2f13
1 parent 5737080
commit 41a2f13
Show file tree

Hide file tree

Showing 19 changed files with 1,008 additions and 0 deletions.
diff --git a/REGIONS/africa.m b/REGIONS/africa.m
@@ -0,0 +1,41 @@
+function varargout=africa(res,buf) 
+% XY=AFRICA(res,buf)
+% AFRICA(...) % Only makes a plot
+%
+% Finds the coordinates of Africa, potentially buffered by some amount.
+%
+% INPUT:
+%
+% res      0 The standard, default values
+%          N Splined values at N times the resolution
+% buf      Distance in degrees that the region outline will be enlarged
+%          by BUFFERM, not necessarily integer, possibly negative
+%          [default: 0]
+%
+% OUTPUT:
+%
+% XY       Closed-curved coordinates of the continent
+%
+% Last modified by charig-at-princeton.edu, 11/23/2011
+% Last modified by fjsimons-at-alum.mit.edu, 11/23/2011
+
+defval('res',0)
+defval('buf',0)
+
+% Parameters that make this the region in question
+regn='africa';
+c11=[0 37.5 ; 342.5 36.5];
+cmn=[52 -35 ; 360 3.5];
+xunt=[261:307 110:260 339:399];
+ofs=[360 0];
+
+% Do it! Make it, load it, save it
+XY=regselect(regn,c11,cmn,xunt,res,buf,ofs);
+
+if nargout==0
+  plot(XY(:,1),XY(:,2),'k-'); axis image; grid on
+end
+
+% Prepare optional output
+varns={XY};
+varargout=varns(1:nargout);
diff --git a/REGIONS/alloceans.m b/REGIONS/alloceans.m
@@ -0,0 +1,115 @@
+function varargout=alloceans(res)
+% XY=ALLOCEANS(res)
+% ALLOCEANS(...) % Only makes a plot
+%
+% Finds the coordinates of some of the worlds' oceans so we can combine
+% its localization kernel with those for the missing continents to turn
+% into Slepian eigenfunctions for all of the world's oceans. Also, there
+% is a polar gap taken out. 
+%
+% INPUT:
+%
+% res      0 The standard, default values
+%          N Splined values at N times the resolution
+% inclang  Inclination angle, determines size of polar caps
+%
+% OUTPUT:
+%
+% XY       Closed-curved coordinates of the world's ocean
+%
+% SEE ALSO: EQBAND
+%
+% Written by D.C. Slobbe, 21/10/2010
+% Last modified by fjsimons-at-alum.mit.edu, 10/08/2011
+
+defval('res',0)
+defval('inclang',66.15);
+
+if res==0
+  fnpl=fullfile(getenv('IFILES'),'COASTS','alloceans.mat');
+else
+  fnpl=fullfile(getenv('IFILES'),'COASTS',sprintf('alloceans-%i.mat',res));
+end
+
+if exist(fnpl,'file')==2 
+  load(fnpl)
+  if nargout==0
+    plot(XY(:,1),XY(:,2),'k-'); axis equal; grid on
+  else
+    varns={XY};
+    varargout=varns(1:nargout);
+  end
+else
+  if res==0
+    ShftAngl    = 180;
+    XYbos       = eqband(0,inclang);
+    XYbos(:,1)  = XYbos(:,1)-360-ShftAngl;
+
+    XY_NA       = namerica(0);
+    XY_NA(:,1)  = XY_NA(:,1)-360;
+
+    XY_GR       = greenland(0);
+    XY_GR(:,1)  = XY_GR(:,1)-360;
+
+    XY_EA       = eurasia(0);
+    XY_EA(:,1)  = XY_EA(:,1)-360;
+
+    % Antarctica receives special treatment
+    clf
+    [~,handl,XYant]=plotcont([0 -62],[360 -83],[],0);
+    delete(handl)
+
+    % Get rid of common NaNs
+    XYant      = XYant(~isnan(XYant(:,1)) & ~isnan(XYant(:,2)),:);
+
+    % Form closed contour
+    XYant      = XYant(~isnan(XYant(:,1)) & ~isnan(XYant(:,2)),:);
+    XYant      = [[0,-90];[0,XYant(227,2)];...
+		  XYant([231:end,1:227],:);[360,-90];[0,-90]];
+    [XY_NA(:,1),XY_NA(:,2)] = poly2cw(XY_NA(:,1), XY_NA(:,2));
+    [XY_GR(:,1),XY_GR(:,2)] = poly2cw(XY_GR(:,1), XY_GR(:,2));
+    [XY_EA(:,1),XY_EA(:,2)] = poly2cw(XY_EA(:,1), XY_EA(:,2));
+    [XYant(:,1),XYant(:,2)] = poly2cw(XYant(:,1), XYant(:,2));
+
+    % Now combine these into an overall interior coastline
+    warning off map:vectorsToGPC:noExternalContours
+    [Xu,Yu]=polybool('subtraction',XYbos(:,1),XYbos(:,2),...
+		     XY_NA(:,1),XY_NA(:,2)); ...
+    [Xu,Yu]=polybool('subtraction',Xu,Yu,XY_GR(:,1),XY_GR(:,2));
+    [Xu,Yu]=polybool('subtraction',Xu,Yu,XYant(:,1),XYant(:,2));
+    [Xu,Yu]=polybool('subtraction',Xu,Yu,XY_EA(:,1),XY_EA(:,2));
+    % Make sure you get both pieces
+    XY_EA(:,1)=XY_EA(:,1)-360;
+    [Xu,Yu]= polybool('subtraction',Xu,Yu,XY_EA(:,1),XY_EA(:,2));
+    XYant(:,1)=XYant(:,1)-360;
+    [Xu,Yu]=polybool('subtraction',Xu,Yu,XYant(:,1),XYant(:,2));
+    XY=[Xu,Yu];
+    warning on map:vectorsToGPC:noExternalContours
+
+    % This should be done manually
+    StrtIndpol  = [0;find(isnan(XY(:,1)));length(XY(:,1))+1];
+    [~,maxdi]   = max(diff(StrtIndpol));
+    XY          = [XY(StrtIndpol(maxdi)+1:StrtIndpol(maxdi+1)-1,:)];
+    [LATo,LONo] = interpm(XY(:,2),XY(:,1),0.25);
+    XY          = [LONo,LATo];
+
+    XY(:,1)     = XY(:,1)+360;
+
+    % Eyeball
+    plot(XY(:,1),XY(:,2),'LineW',2,'Color','k');
+    axis equal
+
+    % Check this out %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    curvecheck(XY(:,1),XY(:,2)) 
+
+    % Check this out %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    save(fnpl,'XY')
+  else
+    XY=alloceans(0);
+    XYb=bezier(XY,res);
+    XY=XYb;
+    save(fnpl,'XY')
+  end
+  varns={XY};
+  varargout=varns(1:nargout);
+end
diff --git a/REGIONS/amazon.m b/REGIONS/amazon.m
@@ -0,0 +1,44 @@
+function varargout=amazon(res,buf)  
+% XY=AMAZON(res,buf)
+% AMAZON(...) % Only makes a plot
+%
+% Finds the coordinates of Amazon, potentially buffered by some amount.
+%
+% INPUT:
+%
+% res      0 The standard, default values
+%          N Splined values at N times the resolution
+% buf      Distance in degrees that the region outline will be enlarged
+%          by BUFFERM, not necessarily integer, possibly negative
+%          [default: 0]
+%
+% OUTPUT:
+%
+% XY       Closed-curved coordinates of the continent
+%
+% Last modified by charig-at-princeton.edu, 11/23/2011
+% Last modified by fjsimons-at-alum.mit.edu, 11/23/2011
+
+defval('res',0)
+defval('buf',0)
+
+% Parameters that make this the region in question
+regn='amazon';
+c11=[280 5];
+cmn=[310 -25];
+xunt=[];
+
+% This admittedly is a special case
+XY=load(fullfile(getenv('IFILES'),'COASTS','amazon'));
+
+% Do it! Make it, load it, save it
+XY=regselect(regn,XY(:,1),XY(:,2),xunt,res,buf);
+
+if nargout==0
+  plot(XY(:,1),XY(:,2),'k-'); axis image; grid on
+end
+
+% Prepare optional output
+varns={XY};
+varargout=varns(1:nargout);
+
diff --git a/REGIONS/antarctica.m b/REGIONS/antarctica.m
@@ -0,0 +1,86 @@
+function varargout=antarctica(res,buf,rotb)
+% [XY,lonc,latc]=antarctica(res,buf,rotb)
+% ANTARCTICA(...) % Only makes a plot
+%
+% Finds the coordinates of Antarctica, but rotates them to an equatorial
+% location so KERNELC doesn't choke on the calculation.
+%
+% INPUT:
+%
+% res      0 The standard, default values
+%          N Splined values at N times the resolution
+%          'rotated' to determine if this region is rotated (it is). We
+%          use this in e.g. GLMALPHA to determine whether to rotate the
+%          eigenfunctions after finding them. You get a logical as output.
+% buf      Distance in degrees by which the region outline will be enlarged
+%          by BUFFERM, not necessarily integer, possibly negative
+%          [default: 0]
+% rotb     0 You want the coordinates on the equator for mathematical 
+%            operations (e.g. you need to make a kernel) [default]
+%          1 You want the coordinates rotated back to their original
+%            location (e.g. for plotting)
+%
+% OUTPUT:
+%
+% XY       Closed-curved coordinates of the continent
+% lonc     The amount by which you need to rotate it back over z
+% latc     The amount by which you need to rotate it back over y
+%
+% See also PLM2ROT, GEOBOXCAP, KLMLMP2ROT, GLMALPHA
+% 
+% Last modified by charig@princeton.edu, 06/24/2016
+% Last modified by fjsimons@princeton.edu, 07/01/2016
+
+defval('res',0)
+defval('buf',0)
+defval('rotb',0)
+
+if ~isstr(res)
+
+  % Some setup
+  c11=[0 -62];
+  cmn=[360 -83];
+
+  % Do it! Make it, load it, save it
+  XY=regselect('antarctica',c11,cmn,[],res,buf);
+
+  % Do some extra here to get the lonc latc
+  [~,~,XY2]=plotcont(c11(:),cmn(:),[],0);
+  close
+  % Get rid of common NaNs
+  XY2=XY2(~isnan(XY2(:,1)) & ~isnan(XY2(:,2)),:);
+  % Get rid of common NaNs
+  XY2=XY2(~isnan(XY2(:,1)) & ~isnan(XY2(:,2)),:);
+  % Find the geographical center and the area to use for rotation
+  % Note: this is based on the unbuffered region and does not change
+  [lonc,latc,A]=rcenter([XY2(:,1) XY2(:,2)]);
+  lonc=-45;
+
+  % Do we return rotated coordinates?
+  if rotb==1
+     [thetap,phip,rotmats]=rottp((90-XY(:,2))*pi/180,XY(:,1)/180*pi,...
+				 -lonc*pi/180,latc*pi/180,0);
+     XY = [phip*180/pi 90-thetap*180/pi];
+  end
+
+  % Make sure the coordinates make sense
+  XY(:,1)=XY(:,1)-360*[XY(:,1)>360];
+  XY(:,1)=XY(:,1)+360*[XY(:,1)<0];
+
+  % Prepare Output
+  varns={XY,lonc,latc};
+  varargout=varns(1:nargout);
+
+elseif strcmp(res,'rotated')
+    % Return a 1 flag as output, indicating the region is a rotated region
+    varargout={1};
+elseif strcmp(res,'demo1')
+    % Make a simple plot
+    XY=antarctica();
+    plot(XY(:,1),XY(:,2),'k-'); axis equal; grid on
+    axis([-30-lonc 30-lonc -30 20])
+end
+
+
+
+
diff --git a/REGIONS/australia.m b/REGIONS/australia.m
@@ -0,0 +1,40 @@
+function varargout=australia(res,buf)  
+% XY=AUSTRALIA(res,buf)
+% AUSTRALIA(...) % Only makes a plot
+%
+% Finds the coordinates of Australia, potentially buffered by some amount.
+%
+% INPUT:
+%
+% res      0 The standard, default values
+%          N Splined values at N times the resolution
+% buf      Distance in degrees that the region outline will be enlarged
+%          by BUFFERM, not necessarily integer, possibly negative
+%          [default: 0]
+%
+% OUTPUT:
+%
+% XY       Closed-curved coordinates of the continent
+%
+% Last modified by charig-at-princeton.edu, 11/23/2011
+% Last modified by fjsimons-at-alum.mit.edu, 11/23/2011
+
+defval('res',0)
+defval('buf',0)
+
+% Parameters that make this the region in question
+regn='australia';
+c11=[112 -10.5];
+cmn=[154 -39];
+xunt=2:166;
+
+% Do it! Make it, load it, save it
+XY=regselect(regn,c11,cmn,xunt,res,buf);
+
+if nargout==0
+  plot(XY(:,1),XY(:,2),'k-'); axis image; grid on
+end
+
+% Prepare optional output
+varns={XY};
+varargout=varns(1:nargout);
diff --git a/REGIONS/baffin.m b/REGIONS/baffin.m
@@ -0,0 +1,39 @@
+function varargout=baffin(res,buf)  
+% XY=baffin(res,buf)
+% baffin(...) % Only makes a plot
+%
+% Finds the coordinates of Baffin Island, potentially buffered by some amount.
+%
+% INPUT:
+%
+% res      0 The standard, default values
+%          N Splined values at N times the resolution
+% buf      Distance in degrees that the region outline will be enlarged
+%          by BUFFERM, not necessarily integer, possibly negative
+%          [default: 0]
+%
+% OUTPUT:
+%
+% XY       Closed-curved coordinates of the continent
+%
+% Last modified by fjsimons-at-alum.mit.edu, 09/23/2014
+
+defval('res',0)
+defval('buf',0)
+
+% Parameters that make this the region in question
+regn='baffin';
+c11=[270 77];
+cmn=[300 58];
+xunt=[14:107];
+
+% Do it! Make it, load it, save it
+XY=regselect(regn,c11,cmn,xunt,res,buf);
+
+if nargout==0
+  plot(XY(:,1),XY(:,2),'k-'); axis equal; grid on
+end
+
+% Prepare optional output
+varns={XY};
+varargout=varns(1:nargout);