![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 26042
- Timestamp:
- 03/05/21 15:37:46 (4 years ago)
- Location:
- issm/trunk-jpl/src/m/contrib/morlighem/modeldata
- Files:
-
- 2 edited
-
interpBedmachineAntarctica.m (modified) (6 diffs)
-
interpBedmachineGreenland.m (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/m/contrib/morlighem/modeldata/interpBedmachineAntarctica.m
r26041 r26042 1 1 function output = interpBedmachineAntarctica(X,Y,string,method,ncdate), 2 %INTERPBEDMACHINEANTARCTICA - interpolate BedMachine data onto X and Y 3 % 4 % Examples: 5 % bed = interpBedmachineAntarctica(X,Y,'bed'); 6 % surface = interpBedmachineAntarctica(X,Y,'surface'); 7 % thickness = interpBedmachineAntarctica(X,Y,'thickness'); 8 % mask = interpBedmachineAntarctica(X,Y,'mask'); 9 % mask = interpBedmachineAntarctica(X,Y,'mask','nearest','../Data/BedMachineAntarctica_2020-07-15_v02.nc'); 10 % 11 % - mask: 0 ocean, 1 land (ice free), 2 grounded ice, 3 floating ice 12 % - source: 1 IBCSO/RTopo-2, 2 MC, 3 interpolation, 4 hydrostatic eq, 13 % 5 Streamline diffusion, 6 Gravity inversion 14 % - optional 4th input argument: interpolation method. 15 % Supported interpolation methos: 'linear','cubic','nearest' 16 % - optional 5th input argument: path to dataset. 17 % 18 % Version 11/30/2018 Mathieu Morlighem mmorligh@uci.edu 2 19 3 20 if nargin<3, string = 'bed'; end 4 21 if nargin<4 5 22 if strcmp(string,'mask') | strcmp(string,'source') 6 7 8 9 23 method='nearest'; % default method 24 else 25 method='cubic'; % default method 26 end 10 27 end 11 28 if nargin<5 12 29 ncdate='2020-07-15'; %BedMachine v2 13 30 end 14 15 31 basename = 'BedMachineAntarctica'; 16 32 17 %read data 18 switch (oshostname()), 19 case {'ronne'} 20 nc=['/home/ModelData/Antarctica/BedMachine/' basename '-' ncdate '.nc']; 21 case {'thwaites','murdo','astrid'} 22 nc=['/u/astrid-r1b/ModelData/BedMachine/' basename '-' ncdate '.nc']; 23 otherwise 24 if nargin==5 25 nc = ncdate; 26 disp(['Using provided path: ' nc]); 27 else 28 error('hostname not supported and path not provided'); 33 if nargin==5 34 ncfile = ncdate; 35 else 36 %List of common paths to try 37 paths = {... 38 ['/u/astrid-r1b/ModelData/BedMachine/' basename '-' ncdate '.nc'],... 39 ['/home/ModelData/Greenland/BedMachine/' basename '-' ncdate '.nc'],... 40 ['/Users/larour/ModelData/BedMachine/' basename '-' ncdate '.nc'],... 41 ['./' basename '-' ncdate '.nc'],... 42 }; 43 44 found = 0; 45 for i=1:numel(paths) 46 if exist(paths{i},'file') 47 ncfile = paths{i}; 48 found = 1; 49 break; 29 50 end 30 end51 end 31 52 32 if nargout==2, 33 string = 'bed'; 53 if ~found 54 error(['Could not find ' basename '-' ncdate '.nc, you can add the path to the list or provide its path as a 5th argument']); 55 end 34 56 end 35 57 36 58 disp([' -- BedMachine Antarctica version: ' ncdate]); 37 xdata = double(ncread(nc ,'x'));38 ydata = double(ncread(nc ,'y'));59 xdata = double(ncread(ncfile,'x')); 60 ydata = double(ncread(ncfile,'y')); 39 61 40 62 offset=2; … … 42 64 xmin=min(X(:)); xmax=max(X(:)); 43 65 posx=find(xdata<=xmax); 66 if isempty(posx), posx=numel(xdata); end 44 67 id1x=max(1,find(xdata>=xmin,1)-offset); 45 68 id2x=min(numel(xdata),posx(end)+offset); … … 47 70 ymin=min(Y(:)); ymax=max(Y(:)); 48 71 posy=find(ydata>=ymin); 72 if isempty(posy), posy=numel(ydata); end 49 73 id1y=max(1,find(ydata<=ymax,1)-offset); 50 74 id2y=min(numel(ydata),posy(end)+offset); … … 52 76 if strcmp(string,'icemask'), 53 77 disp([' -- BedMachine Antarctica: loading ' string]); 54 %data = double(ncread(nc ,'mask'))';55 data = double(ncread(nc ,'mask',[id1x id1y],[id2x-id1x+1 id2y-id1y+1],[1 1]))';78 %data = double(ncread(ncfile,'mask'))'; 79 data = double(ncread(ncfile,'mask',[id1x id1y],[id2x-id1x+1 id2y-id1y+1],[1 1]))'; 56 80 xdata=xdata(id1x:id2x); 57 81 ydata=ydata(id1y:id2y); … … 60 84 else 61 85 disp([' -- BedMachine Antarctica: loading ' string]); 62 %data = double(ncread(nc ,string))';63 data = double(ncread(nc ,string,[id1x id1y],[id2x-id1x+1 id2y-id1y+1],[1 1]))';86 %data = double(ncread(ncfile,string))'; 87 data = double(ncread(ncfile,string,[id1x id1y],[id2x-id1x+1 id2y-id1y+1],[1 1]))'; 64 88 xdata=xdata(id1x:id2x); 65 89 ydata=ydata(id1y:id2y); … … 70 94 if strcmp(string,'mask') | strcmp(string,'source'), 71 95 %Need nearest neighbor to avoid interpolation between 0 and 2 72 tic73 96 output = InterpFromGrid(xdata,ydata,data,double(X),double(Y),'nearest'); 74 toc75 97 %tic 76 98 %output = FastInterp(xdata,ydata,data,X,Y,'nearest'); -
issm/trunk-jpl/src/m/contrib/morlighem/modeldata/interpBedmachineGreenland.m
r25794 r26042 1 1 function output = interpBedmachineGreenland(X,Y,string,ncdate), 2 %INTERPBEDMACHINEGREENLAND - interpolate BedMachine data onto X and Y 3 % 4 % Examples: 5 % bed = interpBedmachineGreenland(X,Y,'bed'); 6 % surface = interpBedmachineGreenland(X,Y,'surface'); 7 % thickness = interpBedmachineGreenland(X,Y,'thickness'); 8 % mask = interpBedmachineGreenland(X,Y,'mask'); 9 % mask = interpBedmachineGreenland(X,Y,'mask','nearest','../Data/BedMachineGreenland_2020-07-15_v03.nc'); 10 % 11 % - mask: 0 ocean, 1 land (ice free), 2 grounded ice, 3 floating ice 12 % - source: 1 IBCSO/RTopo-2, 2 MC, 3 interpolation, 4 hydrostatic eq, 13 % 5 Streamline diffusion, 6 Gravity inversion 14 % - optional 4th input argument: interpolation method. 15 % Supported interpolation methos: 'linear','cubic','nearest' 16 % - optional 5th input argument: path to dataset. 17 % 18 % Version 11/30/2018 Mathieu Morlighem mmorligh@uci.edu 2 19 3 20 if nargin<3, string = 'bed'; end 4 if nargin<4, 21 if nargin<4 22 if strcmp(string,'mask') | strcmp(string,'source') 23 method='nearest'; % default method 24 else 25 method='cubic'; % default method 26 end 27 end 28 if nargin<5 5 29 %ncdate='2015-04-27'; %BedMachine v2 6 30 ncdate='2017-09-25'; %BedMachine v3 7 31 ncdate='2020-04-14'; 8 32 end 9 10 33 basename = 'BedMachineGreenland'; 11 34 12 %List of common paths to try 13 paths = {... 14 ['/u/astrid-r1b/ModelData/ModelData/MCdataset-' ncdate '.nc'],... 15 ['/home/ModelData/Greenland/BedMachine/' basename '-' ncdate '.nc'],... 16 ['/Users/larour/ModelData/BedMachine/' basename '-' ncdate '.nc'],... 17 ['./' basename '-' ncdate '.nc'],... 18 }; 35 if nargin==5 36 ncfile = ncdate; 37 else 38 %List of common paths to try 39 paths = {... 40 ['/u/astrid-r1b/ModelData/ModelData/MCdataset-' ncdate '.nc'],... 41 ['/home/ModelData/Greenland/BedMachine/' basename '-' ncdate '.nc'],... 42 ['/Users/larour/ModelData/BedMachine/' basename '-' ncdate '.nc'],... 43 ['./' basename '-' ncdate '.nc'],... 44 }; 19 45 20 found = 0; 21 for i=1:numel(paths) 22 if exist(paths{i},'file') 23 ncfile = paths{i}; 24 found = 1; 25 break; 46 found = 0; 47 for i=1:numel(paths) 48 if exist(paths{i},'file') 49 ncfile = paths{i}; 50 found = 1; 51 break; 52 end 53 end 54 55 if ~found 56 error(['Could not find ' basename '-' ncdate '.nc, you can add the path to the list or provide its path as a 5th argument']); 26 57 end 27 58 end 28 29 if ~found30 error(['Could not find ' basename '-' ncdate '.nc, you can add the path to the list']);31 end32 33 59 34 60 disp([' -- BedMachine Greenland version: ' ncdate]); … … 57 83 58 84 disp([' -- BedMachine Greenland: interpolating ' string]); 85 disp([' -- Interpolation method: ' method]); 59 86 if strcmp(string,'mask') | strcmp(string,'source'), 60 87 %Need nearest neighbor to avoid interpolation between 0 and 2 … … 64 91 end 65 92 66 %TEST https://www.mathworks.com/matlabcentral/fileexchange/10772-fast-2-dimensional-interpolation 93 end 94 function zi = FastInterp(x,y,data,xi,yi,method) 95 96 %get data size 97 [M N] = size(data); 98 99 % Get X and Y library array spacing 100 ndx = 1/(x(2)-x(1)); ndy = 1/(y(2)-y(1)); 101 % Begin mapping xi and yi vectors onto index space by subtracting library 102 % array minima and scaling to index spacing 103 104 xi = (xi - x(1))*ndx; yi = (yi - y(1))*ndy; 105 106 % Fill Zi with NaNs 107 zi = NaN(size(xi)); 108 109 if strcmpi(method,'nearest'), 110 % Find the nearest point in index space 111 rxi = round(xi)+1; ryi = round(yi)+1; 112 % Find points that are in X,Y range 113 flag = rxi>0 & rxi<=N & ~isnan(rxi) & ryi>0 & ryi<=M & ~isnan(ryi); 114 % Map subscripts to indices 115 ind = ryi + M*(rxi-1); 116 zi(flag) = data(ind(flag)); 117 118 else %Bilinear 119 120 % Transform to unit square 121 fxi = floor(xi)+1; fyi = floor(yi)+1; % x_i and y_i 122 dfxi = xi-fxi+1; dfyi = yi-fyi+1; % Location in unit square 123 124 % flagIn determines whether the requested location is inside of the data arrays 125 flagIn = fxi>0 & fxi<N & ~isnan(fxi) & fyi>0 & fyi<M & ~isnan(fyi); 126 127 %Toss all out-of-bounds variables now to save time 128 fxi = fxi(flagIn); fyi = fyi(flagIn); 129 dfxi = dfxi(flagIn); dfyi = dfyi(flagIn); 130 131 %Find bounding vertices 132 ind1 = fyi + M*(fxi-1); % indices of ( x_i , y_i ) 133 ind2 = fyi + M*fxi; % indices of ( x_i+1 , y_i ) 134 ind3 = fyi + 1 + M*fxi; % indices of ( x_i+1 , y_i+1 ) 135 ind4 = fyi + 1 + M*(fxi-1); % indices of ( x_i , y_i+1 ) 136 137 % Bilinear interpolation 138 zi(flagIn) = ... 139 data(ind1).*(1-dfxi).*(1-dfyi) + ... 140 data(ind2).*dfxi.*(1-dfyi) + ... 141 data(ind4).*(1-dfxi).*dfyi + ... 142 data(ind3).*dfxi.*dfyi; 143 end 144 end
Note:
See TracChangeset
for help on using the changeset viewer.
