Index: /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpBamber2001.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpBamber2001.m (revision 23873) +++ /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpBamber2001.m (revision 23873) @@ -0,0 +1,33 @@ +function [bedout thicknessout] = interpBamber2001(X,Y), + +switch oshostname(), + case {'murdo','thwaites','astrid'} + bamber2001bedpath ='/u/astrid-r1b/ModelData/BamberDEMGreenland5km/bedrock.mat'; + bamber2001thxpath ='/u/astrid-r1b/ModelData/BamberDEMGreenland5km/thickness.mat'; + case {'ronne'} + bamber2001bedpath ='/home/ModelData/Greenland/Bamber2001/bedrock.mat'; + bamber2001thxpath ='/home/ModelData/Greenland/Bamber2001/thickness.mat'; + otherwise + error('machine not supported yet'); +end + +verbose = 0; + +%Convert to Bamber's projections +if verbose, disp(' -- Bamber2001: converting coordinates'); end +[LAT, LON ] = xy2ll(double(X(:)),double(Y(:)),+1,45,70); +[x3971,y3971] = ll2xy(LAT,LON ,+1,39,71); + +if verbose, disp(' -- Bamber2001: loading bed'); end +load(bamber2001bedpath); +if verbose, disp(' -- Bamber2001: interpolating bed'); end +bedout = InterpFromGrid((x_m(1:end-1)+x_m(2:end))/2,(y_m(1:end-1)+y_m(2:end))/2,bedrock,x3971,y3971); +bedout = reshape(bedout,size(X,1),size(X,2)); + +if nargout>1 + if verbose, disp(' -- Bamber2001: loading thickness'); end + load(bamber2001thxpath); + if verbose, disp(' -- Bamber2001: interpolating thickness'); end + thicknessout = InterpFromGrid((x_m(1:end-1)+x_m(2:end))/2,(y_m(1:end-1)+y_m(2:end))/2,thickness,x3971,y3971); + thicknessout = reshape(thicknessout,size(X,1),size(X,2)); +end Index: /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpBamber2013.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpBamber2013.m (revision 23873) +++ /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpBamber2013.m (revision 23873) @@ -0,0 +1,50 @@ +function output = interpBamber2013(X,Y,string), +%INTERPBAMBER2013 - interpolate Bamber 2013 data +% +% Available data: +% BedrockElevation +% SurfaceElevation +% IceThickness +% SurfaceRMSE +% BedrockError +% LandMask (Land mask, 0=ocean, 1=land, 2=ice sheet, 3=non-Greenlandic land, 4=ice shelf) +% NumberAirbornePoints +% Geoid +% BedrockChangeMask +% IceShelfSourceMask +% BedrockElevation_unprocessed +% IceThickness_unprocessed +% BathymetryDataMask + +switch oshostname(), + case {'murdo','thwaites','astrid'} + bamber2013nc='/u/astrid-r1b/morlighe/issmjpl/proj-morlighem/DatasetGreenland/Data/Bamber2013/Greenland_bedrock_topography_V3.nc'; + case {'ronne'} + bamber2013nc='/home/ModelData/Greenland/Bamber2013/Greenland_bedrock_topography_V3.nc'; + otherwise + error('machine not supported yet'); +end +verbose = 0; + +if nargin==2, + string = 'BedrockElevation'; +end + +%Convert to Bamber's projections +if verbose, disp(' -- Bamber2013: converting coordinates'); end +[LAT, LON ] = xy2ll(double(X(:)),double(Y(:)),+1,45,70); +[x3971,y3971] = ll2xy(LAT,LON ,+1,39,71); + +if verbose, disp(' -- Bamber2013: loading coordinates'); end +xdata = double(ncread(bamber2013nc,'projection_x_coordinate'));%*1000; +ydata = double(ncread(bamber2013nc,'projection_y_coordinate'));%*1000; + +if verbose, disp([' -- Bamber2013: loading ' string]); end +data = double(ncread(bamber2013nc,string))'; +if verbose, disp([' -- Bamber2013: interpolating ' string]); end +if strcmpi(string,'LandMask'); + output = InterpFromGrid(xdata,ydata,data,x3971,y3971,'nearest'); +else + output = InterpFromGrid(xdata,ydata,data,x3971,y3971); +end +output = reshape(output,size(X,1),size(X,2)); Index: /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpBedmachineAntarctica.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpBedmachineAntarctica.m (revision 23873) +++ /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpBedmachineAntarctica.m (revision 23873) @@ -0,0 +1,166 @@ +function output = interpBedmachineAntarctica(X,Y,string,method,ncdate), + +if nargin<3, string = 'bed'; end +if nargin<4 + if strcmp(string,'mask') | strcmp(string,'source') + method='nearest'; % default method + else + method='cubic'; % default method + end +end +if nargin<5 + %ncdate='2014-03-12'; + %ncdate='2014-04-28'; + %ncdate='2014-07-31'; + %ncdate='2014-09-23'; + %ncdate='2015-11-06'; + %ncdate='2015-12-09'; + ncdate='2016-06-10'; + ncdate='2016-08-26'; + ncdate='2016-11-17'; + ncdate='2017-01-05'; + ncdate='2017-03-10'; + ncdate='2018-02-02'; + ncdate='2018-03-20'; + ncdate='2018-05-18'; + ncdate='2018-06-08'; + ncdate='2018-09-26'; + ncdate='2018-11-14'; + ncdate='2019-01-24'; + ncdate='2019-01-30'; + ncdate='2019-02-15'; + ncdate='2019-02-21'; % no need of firn correction + ncdate='2019-02-22'; % no need of firn correction + ncdate='2019-02-24'; % no need of firn correction + ncdate='2019-04-15'; % no need of firn correction +end + +date1 = sscanf(ncdate,'%d-%d-%d'); +date2 = datetime(date1(1),date1(2),date1(3)); + +if date2=xmin,1)-offset); +id2x=min(numel(xdata),posx(end)+offset); + +ymin=min(Y(:)); ymax=max(Y(:)); +posy=find(ydata>=ymin); +id1y=max(1,find(ydata<=ymax,1)-offset); +id2y=min(numel(ydata),posy(end)+offset); + +if strcmp(string,'icemask'), + disp([' -- BedMachine Antarctica: loading ' string]); + %data = double(ncread(nc,'mask'))'; + data = double(ncread(nc,'mask',[id1x id1y],[id2x-id1x+1 id2y-id1y+1],[1 1]))'; + xdata=xdata(id1x:id2x); + ydata=ydata(id1y:id2y); + %ice ocean interface is between 0 and 3, so we might get some 1 by interpolating + data(find(data==3))=0; +else + disp([' -- BedMachine Antarctica: loading ' string]); + %data = double(ncread(nc,string))'; + data = double(ncread(nc,string,[id1x id1y],[id2x-id1x+1 id2y-id1y+1],[1 1]))'; + xdata=xdata(id1x:id2x); + ydata=ydata(id1y:id2y); +end + +disp([' -- BedMachine Antarctica: interpolating ' string]); +disp([' -- Interpolation method: ' method]); +if strcmp(string,'mask') | strcmp(string,'source'), + %Need nearest neighbor to avoid interpolation between 0 and 2 + tic + output = InterpFromGrid(xdata,ydata,data,double(X),double(Y),'nearest'); + toc + %tic + %output = FastInterp(xdata,ydata,data,X,Y,'nearest'); + %toc +else + %disp('InterpFromGrid'); + %tic + %output = InterpFromGrid(xdata,ydata,data,double(X),double(Y),'cubic'); + output = InterpFromGrid(xdata,ydata,data,double(X),double(Y),method); % now the interpolation method can be defined by the user + %toc + %disp('FastInterp'); + %tic + %output = FastInterp(xdata,ydata,data,X,Y,'bilinear'); + %toc +end + +end +function zi = FastInterp(x,y,data,xi,yi,method) + + %get data size + [M N] = size(data); + + % Get X and Y library array spacing + ndx = 1/(x(2)-x(1)); ndy = 1/(y(2)-y(1)); + % Begin mapping xi and yi vectors onto index space by subtracting library + % array minima and scaling to index spacing + + xi = (xi - x(1))*ndx; yi = (yi - y(1))*ndy; + + % Fill Zi with NaNs + zi = NaN(size(xi)); + + if strcmpi(method,'nearest'), + % Find the nearest point in index space + rxi = round(xi)+1; ryi = round(yi)+1; + % Find points that are in X,Y range + flag = rxi>0 & rxi<=N & ~isnan(rxi) & ryi>0 & ryi<=M & ~isnan(ryi); + % Map subscripts to indices + ind = ryi + M*(rxi-1); + zi(flag) = data(ind(flag)); + + else %Bilinear + + % Transform to unit square + fxi = floor(xi)+1; fyi = floor(yi)+1; % x_i and y_i + dfxi = xi-fxi+1; dfyi = yi-fyi+1; % Location in unit square + + % flagIn determines whether the requested location is inside of the data arrays + flagIn = fxi>0 & fxi0 & fyi=xmin,1)-offset); +id2x=min(numel(xdata),posx(end)+offset); + +ymin=min(Y(:)); ymax=max(Y(:)); +posy=find(ydata>=ymin); +if isempty(posy), posy=numel(ydata); end +id1y=max(1,find(ydata<=ymax,1)-offset); +id2y=min(numel(ydata),posy(end)+offset); + +disp([' -- BedMachine Greenland: loading ' string]); +data = double(ncread(morlighem2013nc,string,[id1x id1y],[id2x-id1x+1 id2y-id1y+1],[1 1]))'; +xdata=xdata(id1x:id2x); +ydata=ydata(id1y:id2y); +data(find(data==-9999))=NaN; + +disp([' -- BedMachine Greenland: interpolating ' string]); +if strcmp(string,'mask') | strcmp(string,'source'), + %Need nearest neighbor to avoid interpolation between 0 and 2 + output = InterpFromGrid(xdata,ydata,data,double(X),double(Y),'nearest'); +else + output = InterpFromGrid(xdata,ydata,data,double(X),double(Y)); +end + +%TEST https://www.mathworks.com/matlabcentral/fileexchange/10772-fast-2-dimensional-interpolation Index: /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpBedmap.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpBedmap.m (revision 23873) +++ /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpBedmap.m (revision 23873) @@ -0,0 +1,27 @@ +function [dataout] = interpBedmap(X,Y,string), +%INTERPBEDMAP - interpolate bedmap data +% +% Available data: +% 1. bed is bed height +% 2. thickness is ice thickness +% +% Usage: +% [dataout] = interpBedmap(X,Y,string) + +path=[jplsvn() '/proj-morlighem/DatasetAntarctica/Data/BedMap/gridded/']; + +if strcmp(string,'bed'), + path = [path '/bed.mat']; + load(path); + x_m =(x_m(2:end)+x_m(1:end-1))/2.; + y_m =(y_m(2:end)+y_m(1:end-1))/2.; + dataout = InterpFromGrid(x_m,y_m,bed,double(X),double(Y)); +elseif strcmp(string,'thickness') + path = [path '/thickness.mat']; + load(path); + x_m =(x_m(2:end)+x_m(1:end-1))/2.; + y_m =(y_m(2:end)+y_m(1:end-1))/2.; + dataout = InterpFromGrid(x_m,y_m,thickness,double(X),double(Y)); +else + error('not supported'); +end Index: /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpBedmap2.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpBedmap2.m (revision 23873) +++ /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpBedmap2.m (revision 23873) @@ -0,0 +1,95 @@ +function [output] = interpBedmap2(X,Y,string), +%INTERPBEDMAP2 - interpolate bedmap2 data +% +% Available data: +% 1. bed is bed height +% 2. surface is surface height +% 3. thickness is ice thickness +% 4. icemask_grounded_and_shelves is a mask file showing the grounding line and the extent of the floating ice shelves +% 5. rockmask is a mask file showing rock outcrops +% 6. lakemask_vostok is a mask file showing the extent of the lake cavity of Lake Vostok +% 7. bed_uncertainty is the bed uncertainty grid shown in figure 12 of the manuscript +% 8. thickness_uncertainty_5km is the thickness uncertainty grid shown in figure 11 of the manuscript +% 9. coverage is a binary grid showing the distribution of ice thickness data used in the grid of ice thickness +% 10. gl04c_geoid_to_wgs84 is the height conversion values (as floating point) used to convert from WGS84 datum heights to +% g104c geoidal heights (to convert back to WGS84, add this grid) +% +% Usage: +% [dataout] = interpBedmap2(X,Y,string) + +nc = '/home/ModelData/Antarctica/BedMap2/bedmap2_bin/Bedmap2.nc'; +nc = '/Users/mmorligh/temp/Bedmap2.nc'; + +xdata = double(ncread(nc,'x')); +ydata = double(ncread(nc,'y')); + +offset=2; + +xmin=min(X(:)); xmax=max(X(:)); +posx=find(xdata<=xmax); +id1x=max(1,find(xdata>=xmin,1)-offset); +id2x=min(numel(xdata),posx(end)+offset); + +ymin=min(Y(:)); ymax=max(Y(:)); +posy=find(ydata>=ymin); +id1y=max(1,find(ydata<=ymax,1)-offset); +id2y=min(numel(ydata),posy(end)+offset); + +data = double(ncread(nc,string,[id1x id1y],[id2x-id1x+1 id2y-id1y+1],[1 1]))'; +xdata=xdata(id1x:id2x); +ydata=ydata(id1y:id2y); + +if ~strcmp(string,'coverage'), + data(find(data==-9999))=NaN; +end + +if strcmpi(string,'icemask_grounded_and_shelves') | strcmpi(string,'rockmask'), + output = InterpFromGrid(xdata,ydata,data,double(X),double(Y),'nearest'); +else + output = InterpFromGrid(xdata,ydata,data,double(X),double(Y)); % linear interpolation is default +end + +return; +% ================================ OLD =============================================== +%read data +path=['/home/ModelData/Antarctica/BedMap2/bedmap2_bin/']; +if strcmp(string,'gl04c_geoid_to_wgs84'), + filepath = [path '/gl04c_geiod_to_wgs84.flt']; +else + filepath = [path '/bedmap2_' string '.flt']; +end +fid=fopen(filepath,'r','l'); +data=fread(fid,[6667,6667],'float32'); +fclose(fid); + +% define grid +if strcmp(string,'thickness_uncertainty_5km'), + ncols =1361; + nrows =1361; + xll =-3401000; + yll =-3402000; + gridsize =5000; +else + ncols =6667; + nrows =6667; + xll =-3333000; + yll =-3333000; + gridsize =1000; +end +x_m=xll+(0:1:ncols-1)'*gridsize; +y_m=yll+(0:1:nrows-1)'*gridsize; + +%Change default to NaN +if ~strcmp(string,'coverage'), + data(find(data==-9999))=NaN; +end + +%rotate 90 degrees clockwise +data = rot90(data); + +%Interpolate +if strcmpi(string,'icemask_grounded_and_shelves') | strcmpi(string,'rockmask'), + dataout = InterpFromGrid(x_m,y_m,data,double(X),double(Y),'nearest'); +else + dataout = InterpFromGrid(x_m,y_m,data,double(X),double(Y)); +end Index: /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpDhdt.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpDhdt.m (revision 23873) +++ /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpDhdt.m (revision 23873) @@ -0,0 +1,53 @@ +function out = interpDhdt(X,Y), + +switch oshostname(), + case {'ronne'} + dhdtpath='/home/ModelData/Greenland/DHDT/dhdt0306.tif'; + otherwise + error('machine not supported yet'); +end + +%convert coordinates: +[lat lon] = xy2ll(X,Y,+1); +[X Y] = ll2utm(lat,lon,24); + +%Get image info +Tinfo = imfinfo(dhdtpath); +N = Tinfo.Width; +M = Tinfo.Height; +dx = Tinfo.ModelPixelScaleTag(1); +dy = Tinfo.ModelPixelScaleTag(2); +minx = Tinfo.ModelTiepointTag(4); +maxy = Tinfo.ModelTiepointTag(5); + +%Generate vectors +xdata = minx + dx/2 + ((0:N-1).*dx); +ydata = maxy - dy/2 - ((M -1:-1:0).*dy); +ydata = fliplr(ydata); + +%Get pixels we are interested in +offset=2; +xmin=min(X(:)); xmax=max(X(:)); +posx=find(xdata<=xmax); +id1x=max(1,find(xdata>=xmin,1)-offset); +id2x=min(numel(xdata),posx(end)+offset); + +if 0, + ymin=min(Y(:)); ymax=max(Y(:)); + posy=find(ydata<=ymax); + id1y=max(1,find(ydata>=ymin,1)-offset); + id2y=min(numel(ydata),posy(end)+offset); +else + ymin=min(Y(:)); ymax=max(Y(:)); + posy=find(ydata>=ymin); + id1y=max(1,find(ydata<=ymax,1)-offset); + id2y=min(numel(ydata),posy(end)+offset); +end + +data = double(imread(dhdtpath,'PixelRegion',{[id1y,id2y],[id1x,id2x]})); +xdata=xdata(id1x:id2x); +ydata=ydata(id1y:id2y); +data(find(data>+10^3)) = 0; +data(find(data<-10^3)) = 0; + +out = InterpFromGrid(xdata,ydata,data,X,Y); Index: /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpFromGeotiff.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpFromGeotiff.m (revision 23873) +++ /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpFromGeotiff.m (revision 23873) @@ -0,0 +1,40 @@ +function dataout = interpFromGeotiff(geotiffname,X,Y), + +usemap = 0; +if license('test','map_toolbox')==0, + disp('WARNING: map toolbox not installed, trying house code'); + usemap = 0; +elseif license('checkout','map_toolbox')==0 + disp('WARNING: map toolbox not available (checkout failed), trying house code'); + usemap = 0; +end + +if usemap, + [data,R] = geotiffread(geotiffname); + data=double(flipud(data)); + xdata=R.XLimWorld(1):R.DeltaX:R.XLimWorld(2); xdata=xdata(:); + xdata =(xdata(1:end-1)+xdata(2:end))/2; + ydata=R.YLimWorld(2):R.DeltaY:R.YLimWorld(1); ydata=flipud(ydata(:)); + ydata =(ydata(1:end-1)+ydata(2:end))/2; +else + + %Get image info + Tinfo = imfinfo(geotiffname); + N = Tinfo.Width; + M = Tinfo.Height; + dx = Tinfo.ModelPixelScaleTag(1); + dy = Tinfo.ModelPixelScaleTag(2); + minx = Tinfo.ModelTiepointTag(4); + maxy = Tinfo.ModelTiepointTag(5); + + %Generate vectors + xdata = minx + dx/2 + ((0:N-1).*dx); + ydata = maxy - dy/2 - ((M -1:-1:0).*dy); + + %Read image + data=double(flipud(imread(geotiffname))); + data(find(abs(data)>10^30))=NaN; +end + +dataout = InterpFromGrid(xdata,ydata,data,X,Y); +dataout(dataout==-9999)=NaN; Index: /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpGeoid.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpGeoid.m (revision 23873) +++ /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpGeoid.m (revision 23873) @@ -0,0 +1,34 @@ +function [geoid] = interpGeoid(X,Y,varargin), + +switch oshostname(), + case {'murdo','thwaites','astrid'} + rootname=[jplsvn() '/proj-morlighem/DatasetGreenland/Data/Geoid/eigen-6c4-1970.mat']; + case {'ronne'} + rootname='/home/ModelData/Global/Geoid/eigen-6c4-1970.mat'; + otherwise + error('machine not supported yet'); +end +verbose = 1; + +if nargin==3, + hemisphere = varargin{1}; +else + hemisphere = +1; +end + +if hemisphere==+1, + if verbose, disp(' -- Geoid: convert to lat/lon using Greenland projection'); end + [LAT, LON ] = xy2ll(double(X(:)),double(Y(:)),+1,45,70); +else + if verbose, disp(' -- Geoid: convert to lat/lon using Antarctica projection'); end + [LAT, LON ] = xy2ll(double(X(:)),double(Y(:)),-1,0,71); +end +pos=find(LON<0); +LON(pos) =360+LON(pos); + +if verbose, disp(' -- Geoid: loading eigen-6c4 '); end +A=load(rootname); + +if verbose, disp(' -- Geoid: interpolating'); end +geoid = InterpFromGrid(A.lon,A.lat,A.geoid,LON,LAT); +geoid = reshape(geoid,size(X)); Index: /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpGimpdem.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpGimpdem.m (revision 23873) +++ /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpGimpdem.m (revision 23873) @@ -0,0 +1,68 @@ +function sout = interpGimpdem(X,Y), + +switch oshostname(), + case {'murdo','thwaites','astrid'} + howatpath='/u/astrid-r1b/morlighe/issmjpl/proj-morlighem/DatasetGreenland/Data/gimpdem/gimpdem_90m.tif'; + case {'ronne'} + howatpath='/home/ModelData/Greenland/gimpdem/gimpdem_90m.tif'; + otherwise + error('machine not supported yet'); +end + +usemap = 0; +if license('test','map_toolbox')==0, + disp('WARNING: map toolbox not installed, trying house code'); + usemap = 0; +elseif license('checkout','map_toolbox')==0 + disp('WARNING: map toolbox not available (checkout failed), trying house code'); + usemap = 0; +end + +if usemap, + [data,R] = geotiffread(howatpath); + data=double(flipud(data)); + xdata=R.XLimWorld(1):R.DeltaX:R.XLimWorld(2); xdata=xdata(:); + xdata =(xdata(1:end-1)+xdata(2:end))/2; + ydata=R.YLimWorld(2):R.DeltaY:R.YLimWorld(1); ydata=flipud(ydata(:)); + ydata =(ydata(1:end-1)+ydata(2:end))/2; +else + + %Get image info + Tinfo = imfinfo(howatpath); + N = Tinfo.Width; + M = Tinfo.Height; + dx = Tinfo.ModelPixelScaleTag(1); + dy = Tinfo.ModelPixelScaleTag(2); + minx = Tinfo.ModelTiepointTag(4); + maxy = Tinfo.ModelTiepointTag(5); + + %Generate vectors + xdata = minx + dx/2 + ((0:N-1).*dx); + ydata = maxy - dy/2 - ((M -1:-1:0).*dy); + ydata = fliplr(ydata); + + %Get pixels we are interested in + offset=2; + xmin=min(X(:)); xmax=max(X(:)); + posx=find(xdata<=xmax); + id1x=max(1,find(xdata>=xmin,1)-offset); + id2x=min(numel(xdata),posx(end)+offset); + + if 0, + ymin=min(Y(:)); ymax=max(Y(:)); + posy=find(ydata<=ymax); + id1y=max(1,find(ydata>=ymin,1)-offset); + id2y=min(numel(ydata),posy(end)+offset); + else + ymin=min(Y(:)); ymax=max(Y(:)); + posy=find(ydata>=ymin); + id1y=max(1,find(ydata<=ymax,1)-offset); + id2y=min(numel(ydata),posy(end)+offset); + end + + data = double(imread(howatpath,'PixelRegion',{[id1y,id2y],[id1x,id2x]})); + xdata=xdata(id1x:id2x); + ydata=ydata(id1y:id2y); +end + +sout = InterpFromGrid(xdata,ydata,data,X,Y); Index: /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpGimpicemask.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpGimpicemask.m (revision 23873) +++ /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpGimpicemask.m (revision 23873) @@ -0,0 +1,69 @@ +function sout = interpGimpicemask(X,Y), + +switch oshostname(), + case {'ronne'} + howatpath='/home/ModelData/Greenland/gimpmask/GimpIceMask_90m.tif'; + otherwise + error('machine not supported yet'); +end + +usemap = 0; +if license('test','map_toolbox')==0, + disp('WARNING: map toolbox not installed, trying house code'); + usemap = 0; +elseif license('checkout','map_toolbox')==0 + disp('WARNING: map toolbox not available (checkout failed), trying house code'); + usemap = 0; +end + +if usemap, + [data,R] = geotiffread(howatpath); + data=double(flipud(data)); + xdata=R.XLimWorld(1):R.DeltaX:R.XLimWorld(2); xdata=xdata(:); + xdata =(xdata(1:end-1)+xdata(2:end))/2; + ydata=R.YLimWorld(2):R.DeltaY:R.YLimWorld(1); ydata=flipud(ydata(:)); + ydata =(ydata(1:end-1)+ydata(2:end))/2; +else + + %Get image info + Tinfo = imfinfo(howatpath); + N = Tinfo.Width; + M = Tinfo.Height; + dx = Tinfo.ModelPixelScaleTag(1); + dy = Tinfo.ModelPixelScaleTag(2); + minx = Tinfo.ModelTiepointTag(4); + maxy = Tinfo.ModelTiepointTag(5); + + %Generate vectors + xdata = minx + dx/2 + ((0:N-1).*dx); + ydata = maxy - dy/2 - ((M -1:-1:0).*dy); + ydata = fliplr(ydata); + + %Get pixels we are interested in + offset=2; + xmin=min(X(:)); xmax=max(X(:)); + posx=find(xdata<=xmax); + id1x=max(1,find(xdata>=xmin,1)-offset); + id2x=min(numel(xdata),posx(end)+offset); + + if 0, + ymin=min(Y(:)); ymax=max(Y(:)); + posy=find(ydata<=ymax); + id1y=max(1,find(ydata>=ymin,1)-offset); + id2y=min(numel(ydata),posy(end)+offset); + else + ymin=min(Y(:)); ymax=max(Y(:)); + posy=find(ydata>=ymin); + id1y=max(1,find(ydata<=ymax,1)-offset); + id2y=min(numel(ydata),posy(end)+offset); + end + + data = double(imread(howatpath,'PixelRegion',{[id1y,id2y],[id1x,id2x]})); + xdata=xdata(id1x:id2x); + ydata=ydata(id1y:id2y); +end + +sout = InterpFromGrid(xdata,ydata,data,X,Y,'nearest'); + +%Post process output (undefined = not ice) +sout(find(sout==-9999))=0; Index: /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpGimpoceanmask.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpGimpoceanmask.m (revision 23873) +++ /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpGimpoceanmask.m (revision 23873) @@ -0,0 +1,69 @@ +function sout = interpGimpoceanmask(X,Y), + +switch oshostname(), + case {'ronne'} + howatpath='/home/ModelData/Greenland/gimpmask/GimpOceanMask_90m.tif'; + otherwise + error('machine not supported yet'); +end + +usemap = 0; +if license('test','map_toolbox')==0, + disp('WARNING: map toolbox not installed, trying house code'); + usemap = 0; +elseif license('checkout','map_toolbox')==0 + disp('WARNING: map toolbox not available (checkout failed), trying house code'); + usemap = 0; +end + +if usemap, + [data,R] = geotiffread(howatpath); + data=double(flipud(data)); + xdata=R.XLimWorld(1):R.DeltaX:R.XLimWorld(2); xdata=xdata(:); + xdata =(xdata(1:end-1)+xdata(2:end))/2; + ydata=R.YLimWorld(2):R.DeltaY:R.YLimWorld(1); ydata=flipud(ydata(:)); + ydata =(ydata(1:end-1)+ydata(2:end))/2; +else + + %Get image info + Tinfo = imfinfo(howatpath); + N = Tinfo.Width; + M = Tinfo.Height; + dx = Tinfo.ModelPixelScaleTag(1); + dy = Tinfo.ModelPixelScaleTag(2); + minx = Tinfo.ModelTiepointTag(4); + maxy = Tinfo.ModelTiepointTag(5); + + %Generate vectors + xdata = minx + dx/2 + ((0:N-1).*dx); + ydata = maxy - dy/2 - ((M -1:-1:0).*dy); + ydata = fliplr(ydata); + + %Get pixels we are interested in + offset=2; + xmin=min(X(:)); xmax=max(X(:)); + posx=find(xdata<=xmax); + id1x=max(1,find(xdata>=xmin,1)-offset); + id2x=min(numel(xdata),posx(end)+offset); + + if 0, + ymin=min(Y(:)); ymax=max(Y(:)); + posy=find(ydata<=ymax); + id1y=max(1,find(ydata>=ymin,1)-offset); + id2y=min(numel(ydata),posy(end)+offset); + else + ymin=min(Y(:)); ymax=max(Y(:)); + posy=find(ydata>=ymin); + id1y=max(1,find(ydata<=ymax,1)-offset); + id2y=min(numel(ydata),posy(end)+offset); + end + + data = double(imread(howatpath,'PixelRegion',{[id1y,id2y],[id1x,id2x]})); + xdata=xdata(id1x:id2x); + ydata=ydata(id1y:id2y); +end + +sout = InterpFromGrid(xdata,ydata,data,X,Y,'nearest'); + +%Post process output (undefined = ocean) +sout(find(sout==-9999))=1; Index: /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpGridsCReSIS.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpGridsCReSIS.m (revision 23873) +++ /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpGridsCReSIS.m (revision 23873) @@ -0,0 +1,30 @@ +function output = interpGridsCReSIS(X,Y,filename), + +%Convert to lat/lon +disp(' -- Griggs2013: converting coordinates'); +[LAT, LON ] = xy2ll(double(X(:)),double(Y(:)),+1,45,70); + +disp([' -- GridsCReSIS: loading data']); +if ~exist(filename) + error([filename ' does not exist']); +end +fid = fopen(filename); +for i=1:6, + thisline = fgetl(fid); + dummy = regexp(thisline,'(\S+)','match'); + if strcmp(dummy{1},'ncols'), ncols=str2num(dummy{2}); end + if strcmp(dummy{1},'nrows'), nrows=str2num(dummy{2}); end + if strcmp(dummy{1},'xllcorner'), xllcorner=str2num(dummy{2}); end + if strcmp(dummy{1},'yllcorner'), yllcorner=str2num(dummy{2}); end + if strcmp(dummy{1},'cellsize'), cellsize=str2num(dummy{2}); end + if strcmp(dummy{1},'NODATA_value'), nodata=str2num(dummy{2}); end +end +data = fscanf(fid,'%g %g %g %g %g',[ncols nrows])'; +fclose(fid); + +xdata=linspace(xllcorner+cellsize/2,xllcorner+cellsize/2+(ncols-1)*cellsize,ncols); +ydata=linspace(yllcorner+cellsize/2,yllcorner+cellsize/2+(nrows-1)*cellsize,nrows); + +disp([' -- GridsCReSIS: interpolating ']); +output = InterpFromGrid(xdata,ydata,data,LAT,LON); +output = reshape(output,size(X,1),size(X,2)); Index: /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpGriggs2013.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpGriggs2013.m (revision 23873) +++ /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpGriggs2013.m (revision 23873) @@ -0,0 +1,29 @@ +function output = interpGriggs2013(X,Y,string), + +disp('============================================'); +disp(' '); +disp('WARNING: interpBamber2013 should now be used'); +disp(' '); +disp('============================================'); +error('interpBamber2013 should now be used'); +griggs2013nc='/u/astrid-r1b/morlighe/issmjpl/proj-morlighem/DatasetGreenland/Data/Griggs2012/Greenland_bedrock_topography_and_geometry_062012_JGriggs.nc'; +verbose = 0; + +if nargout==2, + string = 'BedrockElevation'; +end + +%Convert to Bamber's projections +if verbose, disp(' -- Griggs2013: converting coordinates'); end +[LAT, LON ] = xy2ll(double(X(:)),double(Y(:)),+1,45,70); +[x3971,y3971] = ll2xy(LAT,LON ,+1,39,71); + +if verbose, disp(' -- Griggs2013: loading coordinates'); end +xdata = double(ncread(griggs2013nc,'projection_x_coordinate'))*1000; +ydata = double(ncread(griggs2013nc,'projection_y_coordinate'))*1000; + +if verbose, disp([' -- Griggs2013: loading ' string]); end +data = double(ncread(griggs2013nc,string))'; +if verbose, disp([' -- Griggs2013: interpolating ' string]); end +output = InterpFromGrid(xdata,ydata,data,x3971,y3971); +output = reshape(output,size(X,1),size(X,2)); Index: /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpIBCSO.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpIBCSO.m (revision 23873) +++ /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpIBCSO.m (revision 23873) @@ -0,0 +1,38 @@ +function [bedout sid] = interpIBCSO(X,Y), + +%read data +switch (oshostname()), + case {'ronne'} + ncpath='/home/ModelData/Antarctica/IBCSO/ibcso_v1_bed.grd'; + sidpath='/home/ModelData/Antarctica/IBCSO/ibcso_v1_sid.grd'; + otherwise + error('hostname not supported yet'); +end + +disp(' -- IBCSO: loading bathymetry'); +x_range = double(ncread(ncpath,'x_range')); +y_range = double(ncread(ncpath,'y_range')); +spacing = double(ncread(ncpath,'spacing')); +xdata = (x_range(1)-spacing(1)/2) : spacing(1) : (x_range(2)-spacing(1)/2); +ydata = (y_range(1)-spacing(2)/2) : spacing(2) : (y_range(2)-spacing(2)/2); +data = double(ncread(ncpath,'z')); +data(find(data==-9999 | isinf(data))) = NaN; +data = reshape(data,[numel(xdata) numel(ydata)])'; +disp(' -- IBCSO: interpolating bed'); +bedout = InterpFromGrid(xdata,fliplr(ydata),data,double(X),double(Y)); + +if nargout==2, + disp(' -- IBCSO: bathymetry sid'); + xdata = ncread(sidpath,'x'); + ydata = ncread(sidpath,'y'); + data = ncread(sidpath,'z')'; + disp(' -- IBCSO: transforming coordinates'); + [LAT,LON] = xy2ll(double(X(:)),double(Y(:)),-1,0,71); + [x065,y065] = ll2xy(LAT,LON,-1,0,65); + x065 = reshape(x065,size(X)); + y065 = reshape(y065,size(Y)); + disp(' -- IBCSO: interpolating sids'); + sid = InterpFromGrid(xdata,ydata,data,x065,y065,'nearest'); + sid(find(sid<200000)) = 0; + sid(find(sid>399999)) = 0; +end Index: /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpJakobsson2012.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpJakobsson2012.m (revision 23873) +++ /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpJakobsson2012.m (revision 23873) @@ -0,0 +1,35 @@ +function [bedout sourceout] = interpJakobsson2012(X,Y,string), + +switch oshostname(), + case {'murdo','thwaites','astrid'} + ncpath ='/u/astrid-r1b/morlighe/issmjpl/proj-morlighem/DatasetGreenland/Data/IBCAO/IBCAO_V3_500m_RR.grd'; + case {'ronne'} + ncpath ='/home/ModelData/Greenland/IBCAO/IBCAO_V3_500m_RR.grd'; + otherwise + error('machine not supported yet'); +end + +%Convert to IBCAO projections +disp(' -- Jakobsson2012: converting coordinates'); +[LAT, LON ] = xy2ll(double(X(:)),double(Y(:)),+1,45,70); +[x0075,y0075] = ll2xy(LAT,LON,+1,0,75); + +disp(' -- Jakobsson2012: loading bathymetry'); +xdata = double(ncread(ncpath,'x')); +ydata = double(ncread(ncpath,'y')); +data = double(ncread(ncpath,'z'))'; + +disp(' -- Jakobsson2012: interpolating bed'); +bedout = InterpFromGrid(xdata,ydata,data,x0075,y0075); +bedout = reshape(bedout,size(X,1),size(X,2)); + +if nargout==2, + ncpath ='/home/ModelData/Greenland/IBCAO/IBCAO_V3_SID_500m.grd'; + disp(' -- Jakobsson2012: loading source'); + xdata = double(ncread(ncpath,'x')); + ydata = double(ncread(ncpath,'y')); + data = double(ncread(ncpath,'z'))'; + disp(' -- Jakobsson2012: interpolating source'); + sourceout = InterpFromGrid(xdata,ydata,data,x0075,y0075,'nearest'); + sourceout = reshape(sourceout,size(X,1),size(X,2)); +end Index: /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpJoughin.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpJoughin.m (revision 23873) +++ /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpJoughin.m (revision 23873) @@ -0,0 +1,97 @@ +function [vxout vyout] = interpJoughin(X,Y,Date), + %Available dates: + % 2000 2005 2006 2007 2008 + +switch oshostname(), + case {'murdo','thwaites','astrid'} + if nargin==3, + rootname = ['/u/astrid-r1b/morlighe/issmjpl/proj-morlighem/DatasetGreenland/Data/Vel/Joughin/' num2str(Date) '/']; + else + error('not supported'); + end + case {'ronne'} + error('not supported'); + otherwise + error('machine not supported yet'); +end +verbose = 1; + +if ~exist(rootname,'dir'), + error(['file ' rootname ' not found']); +end + +rootname = [rootname 'greenland_vel_mosaic500_' num2str(Date) '_' num2str(Date+1)]; + +if verbose, disp(' -- Joughin: loading vx'); end +[data,R] = geotiffread([rootname '_vx.tif']); +pos=find(data<-10^9); data(pos)=NaN; +data=double(flipud(data)); +xdata=R.XLimWorld(1):R.DeltaX:R.XLimWorld(2); xdata=xdata(:); +xdata =(xdata(1:end-1)+xdata(2:end))/2; +ydata=R.YLimWorld(2):R.DeltaY:R.YLimWorld(1); ydata=flipud(ydata(:)); +ydata =(ydata(1:end-1)+ydata(2:end))/2; +if verbose, disp(' -- Joughin: interpolating vx'); end +vxout = InterpFromGrid(xdata,ydata,data,X,Y); +vxout = reshape(vxout,size(X,1),size(X,2)); + +if verbose, disp(' -- Joughin: loading vy'); end +[data,R] = geotiffread([rootname '_vy.tif']); +pos=find(data<-10^9); data(pos)=NaN; +data=double(flipud(data)); +xdata=R.XLimWorld(1):R.DeltaX:R.XLimWorld(2); xdata=xdata(:); +xdata =(xdata(1:end-1)+xdata(2:end))/2; +ydata=R.YLimWorld(2):R.DeltaY:R.YLimWorld(1); ydata=flipud(ydata(:)); +ydata =(ydata(1:end-1)+ydata(2:end))/2; +if verbose, disp(' -- Joughin: interpolating vy'); end +vyout = InterpFromGrid(xdata,ydata,data,X,Y); +vyout = reshape(vyout,size(X,1),size(X,2)); +return + +% Get geodat info +if verbose, disp(' -- Joughin: loading geodat info'); end +xd=readgeodat(strcat(rootname,'.vx.geodat')); +xmin=xd(3,1)*1000.+xd(2,1)/2; +xmax=xd(3,1)*1000.+(xd(2,1)-1)*xd(1,1)+xd(2,1)/2; +ymin=xd(3,2)*1000.+xd(2,2)/2; +ymax=xd(3,2)*1000.+(xd(2,2)-1)*xd(1,2)+xd(2,2)/2; +%xmin=xd(3,1)*1000.; +%xmax=xd(3,1)*1000.+(xd(2,1)-1)*xd(1,1); +%ymin=xd(3,2)*1000.; +%ymax=xd(3,2)*1000.+(xd(2,2)-1)*xd(1,2); +xdata=linspace(xmin,xmax,xd(1,1)); +ydata=linspace(ymin,ymax,xd(1,2)); + +% Vx component +if verbose, disp(' -- Joughin: loading vx'); end +fid = fopen(strcat(rootname,'.vx'),'r','ieee-be'); +[data,count]=fread(fid,[xd(1,1) xd(1,2)],'float32'); +fclose(fid); + +if verbose, disp(' -- Joughin: interpolating vx'); end +vxout = InterpFromGrid(xdata,ydata,data',X,Y); +vxout = reshape(vxout,size(X,1),size(X,2)); + +% Vy component +fid = fopen(strcat(rootname,'.vy'),'r','ieee-be'); +[data,count]=fread(fid,[xd(1,1) xd(1,2)],'float32'); +fclose(fid); +vyout = InterpFromGrid(xdata,ydata,data',X,Y); +vyout = reshape(vyout,size(X,1),size(X,2)); + +end + +function xgeo=readgeodat(filein) +% Read a geodat file +fid = fopen(filein,'r'); +xgeo=zeros(3,2); +i=1; +while ~feof(fid), + line=fgets(fid); + [A,count]=sscanf(line,'%f %f',[1 2]); + if(count == 2) + xgeo(i,:)=A; + i=i+1; + end +end +fclose(fid); +end Index: /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpJoughinCompositeGreenland.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpJoughinCompositeGreenland.m (revision 23873) +++ /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpJoughinCompositeGreenland.m (revision 23873) @@ -0,0 +1,32 @@ +function [vxout vyout] = interpJoughinCompositeGreenland(X,Y), + +%data=load(['/u/astrid-r1b/morlighe/issmjpl/proj-morlighem/DatasetGreenland/Data/VelJoughin/IanGreenVel.mat']); +filename = '/home/ModelData/Greenland/VelJoughin/IanGreenVel.mat'; + +%Figure out what subset of the matrix should be read +load(filename,'x_m','y_m'); +velfile = matfile(filename); + +offset=2; + +xmin=min(X(:)); xmax=max(X(:)); +posx=find(x_m<=xmax); +id1x=max(1,find(x_m>=xmin,1)-offset); +id2x=min(numel(x_m),posx(end)+offset); + +ymin=min(Y(:)); ymax=max(Y(:)); +posy=find(y_m>=ymin); +id1y=max(1,find(y_m<=ymax,1)-offset); +id2y=min(numel(y_m),posy(end)+offset); + +vx = velfile.vx(id1y:id2y,id1x:id2x); +vy = velfile.vy(id1y:id2y,id1x:id2x); +x = x_m(id1x:id2x); +y = y_m(id1y:id2y); + +vxout = InterpFromGrid(x,y,double(vx),X,Y); +vyout = InterpFromGrid(x,y,double(vy),X,Y); + +if nargout==1, + vxout = sqrt(vxout.^2+vyout.^2); +end Index: /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpJoughinMosaic.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpJoughinMosaic.m (revision 23873) +++ /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpJoughinMosaic.m (revision 23873) @@ -0,0 +1,34 @@ +function [vxout vyout] = interpJoughinMosaic(X,Y), + +switch oshostname(), + case {'ronne'} + filename = '/home/ModelData/Greenland/VelJoughin/IanGreenVel.mat'; + otherwise + error('machine not supported yet'); +end +verbose = 1; + +%Figure out what subset of the matrix should be read +load(filename,'x_m','y_m'); +velfile = matfile(filename); + +offset=2; + +xmin=min(X(:)); xmax=max(X(:)); +posx=find(x_m<=xmax); +id1x=max(1,find(x_m>=xmin,1)-offset); +id2x=min(numel(x_m),posx(end)+offset); + +ymin=min(Y(:)); ymax=max(Y(:)); +posy=find(y_m>=ymin); +id1y=max(1,find(y_m<=ymax,1)-offset); +id2y=min(numel(y_m),posy(end)+offset); + +vx = velfile.vx(id1y:id2y,id1x:id2x); +vy = velfile.vy(id1y:id2y,id1x:id2x); +x_m = x_m(id1x:id2x); +y_m = y_m(id1y:id2y); + +%load(filename); +vxout = InterpFromGrid(x_m,y_m,vx,X,Y); +vyout = InterpFromGrid(x_m,y_m,vy,X,Y); Index: /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpMouginotAnt2016.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpMouginotAnt2016.m (revision 23873) +++ /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpMouginotAnt2016.m (revision 23873) @@ -0,0 +1,46 @@ +function [vxout vyout]= interpMouginotAnt2016(X,Y), + +%read data +switch (oshostname()), + case {'ronne'} + filename = '/home/ModelData/Antarctica/MouginotVel/vel_ant_5Apr2016.mat'; + case {'thwaites','murdo','astrid'} + filename = '/u/astrid-r1b/ModelData/RignotAntarcticaVelMosaic450m/vel_ant_5Apr2016.mat'; + otherwise + error('hostname not supported yet'); +end + +%Figure out what subset of the matrix should be read +load(filename,'x','y'); +velfile = matfile(filename); + +offset=2; + +xmin=min(X(:)); xmax=max(X(:)); +posx=find(x<=xmax); +id1x=max(1,find(x>=xmin,1)-offset); +id2x=min(numel(x),posx(end)+offset); + +if y(2)-y(1)<0 + ymin=min(Y(:)); ymax=max(Y(:)); + posy=find(y>=ymin); + id1y=max(1,find(y<=ymax,1)-offset); + id2y=min(numel(y),posy(end)+offset); +else + ymin=min(X(:)); ymax=max(X(:)); + posy=find(y<=ymax); + id1y=max(1,find(y>=ymin,1)-offset); + id2y=min(numel(y),posy(end)+offset); +end + +vx = velfile.vx(id1y:id2y,id1x:id2x); +vy = velfile.vy(id1y:id2y,id1x:id2x); +x = x(id1x:id2x); +y = y(id1y:id2y); + +vxout = InterpFromGrid(x,y,double(vx),X,Y); +vyout = InterpFromGrid(x,y,double(vy),X,Y); + +if nargout==1, + vxout = sqrt(vxout.^2+vyout.^2); +end Index: /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpMouginotAnt2017.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpMouginotAnt2017.m (revision 23873) +++ /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpMouginotAnt2017.m (revision 23873) @@ -0,0 +1,36 @@ +function [vxout vyout]= interpRignot2012(X,Y), + +filename = '/home/ModelData/Greenland/VelMouginot/RignotGreenland2012Vel.mat'; + + +%Figure out what subset of the matrix should be read +load(filename,'x','y'); +velfile = matfile(filename); + +offset=2; + +xmin=min(X(:)); xmax=max(X(:)); +posx=find(x<=xmax); +id1x=max(1,find(x>=xmin,1)-offset); +id2x=min(numel(x),posx(end)+offset); + +ymin=min(Y(:)); ymax=max(Y(:)); +%posy=find(y>=ymin); +%id1y=max(1,find(y<=ymax,1)-offset); +%id2y=min(numel(y),posy(end)+offset); +posy=find(y<=ymax); +id1y=max(1,find(y>=ymin,1)-offset); +id2y=min(numel(y),posy(end)+offset); + +vx = velfile.vx(id1y:id2y,id1x:id2x); +vy = velfile.vy(id1y:id2y,id1x:id2x); +x = x(id1x:id2x); +y = y(id1y:id2y); + +%load(filename); +vxout = InterpFromGrid(x,y,double(vx),X,Y); +vyout = InterpFromGrid(x,y,double(vy),X,Y); + +if nargout==1, + vxout = sqrt(vxout.^2+vyout.^2); +end Index: /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpMouginotAntTimeSeries1973to2018.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpMouginotAntTimeSeries1973to2018.m (revision 23873) +++ /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpMouginotAntTimeSeries1973to2018.m (revision 23873) @@ -0,0 +1,160 @@ +function [vxout vyout]= interpMouginotAntTimeSeries1973to2018(X,Y,T) +%INTERPMOUGINOTANTTIMESERIES1973TO2018 - interpolate observed (time series) velocities +% +% Inputs +% X,Y: spatial (scatter) coordinates +% T: time (indexed by YEAR1 (YEAR2 is optional); see below) +% +% Outputs +% vxout,vyout: interpolated velocities at X,Y, for each time requested in T +% +% Available time series: +% +% YEAR1 YEAR2 +% 1 1973 1975 +% 2 1973 1984 +% 3 1973 1988 +% 4 1984 1988 +% 5 1986 1988 +% 6 1988 1990 +% 7 1991 1992 +% 8 1995 1996 +% 9 2000 2001 +% 10 2002 2003 +% 11 2003 2004 +% 12 2005 2006 +% 13 2006 2007 +% 14 2007 2008 +% 15 2008 2009 +% 16 2009 2010 +% 17 2010 2011 +% 18 2011 2012 +% 19 2012 2013 +% 20 2013 2014 +% 21 2014 2015 +% 22 2015 2016 +% 23 2016 2017 +% 24 2017 2018 +% +% Usage: +% T refers to YEAR1, but the user can also use YEAR2 (e.g., the "1973" case in YEAR1). +% +% Then, these codes generate the same results: +% +% [vxout vyout]= interpMouginotAntTimeSeries1973to2018(md.mesh.x,md.mesh.y,[1986; 1991; 1995; 2000]); +% [vxout vyout]= interpMouginotAntTimeSeries1973to2018(md.mesh.x,md.mesh.y,[1986 1988; 1991 1992; 1995 1996; 2000 2001]); +% +% Another example: +% [vxout vyout]= interpMouginotAntTimeSeries1973to2018(md.mesh.x,md.mesh.y,[1973 1975; 1973 1988; 1991 1992; 2011 2012]); + +%read data +switch (oshostname()), + case {'ronne'} + nc = '/home/ModelData/Antarctica/MouginotVel/ASE_TimeSeries_1973-2018.nc'; + otherwise + error('hostname not supported yet'); +end + +xdata = double(ncread(nc,'x')); +ydata = double(ncread(nc,'y')); +year1 = ncread(nc,'YEAR1'); +year2 = ncread(nc,'YEAR2'); + +% get the positions related to T +if nargin==3 + % initial checks %{{{ + if size(T,2)>2 | size(T,1)<1 | size(T,2)<1, + error('Size of input T not supported!'); + end + if size(T,2)==1 & any(T(:,1)==1973), + disp(' '); + disp(' Found year=1973 in T (array). Please, specify the data series using a second index.'); + disp(' Data available for 1973:'); + disp(' 1973 1975'); + disp(' 1973 1984'); + disp(' 1973 1988'); + disp(' '); + disp(' Usage:'); + disp(' [vxout vyout]= interpMouginotAntTimeSeries1973to2018(md.mesh.x,md.mesh.y,[1973 1975; 1973 1988; 1991 1992; 2011 2012])'); + disp(' '); + error(' Change input T before continuing.'); + end %}}} + pos = []; + for i=1:size(T,1), + flag = (T(i,1)==year1); + if size(T,2)==2, % ok, check both indexes (year1 and year2) + flag = (T(i,1)==year1).*(T(i,2)==year2); + end + pos = [pos; find(flag)]; + end + % check again {{{ + if length(pos)~=size(T,1) | length(unique(pos))~=length(pos), + disp(' '); + disp(' Time resquested does not exist in data set or is repeated!'); + disp(' Data resquested:'); + for i=1:length(T(:,1)), + str = [' ' int2str(T(i,1)) ' ']; + if size(T,2)==2, % ok, check both indexes (year1 and year2) + str = [str int2str(T(i,2))]; + end + disp(str); + end + disp(' '); + disp(' Data available (24 series):'); + for i=1:length(year1), + str = [' ' int2str(year1(i)) ' ' int2str(year2(i))]; + disp(str); + end + disp(' '); + disp(' Usage:'); + disp(' [vxout vyout]= interpMouginotAntTimeSeries1973to2018(md.mesh.x,md.mesh.y,[1986; 1991; 1995; 2000])'); + disp(' [vxout vyout]= interpMouginotAntTimeSeries1973to2018(md.mesh.x,md.mesh.y,[1986 1988; 1991 1992; 1995 1996; 2000 2001])'); + disp(' [vxout vyout]= interpMouginotAntTimeSeries1973to2018(md.mesh.x,md.mesh.y,[1973 1975; 1973 1988; 1991 1992; 2011 2012])'); + disp(' '); + error(' Change input T before continuing.'); + end%}}} +elseif nargin<3, + pos = 1:24; % all available data +else + error('nargin not supported yet!'); +end + +% get the spatial positions +offset=2; + +xmin=min(X(:)); xmax=max(X(:)); +posx=find(xdata<=xmax); +id1x=max(1,find(xdata>=xmin,1)-offset); +id2x=min(numel(xdata),posx(end)+offset); + +ymin=min(Y(:)); ymax=max(Y(:)); +posy=find(ydata>=ymin); +id1y=max(1,find(ydata<=ymax,1)-offset); +id2y=min(numel(ydata),posy(end)+offset); + +disp([' -- Mouginot Time Series 1973 to 2018: loading velocities']); +vxdata = []; +vydata = []; +for i=1:length(pos), + disp([' step = ' int2str(i) '/' int2str(length(pos)) ', position = ' int2str(pos(i)) ', year = ' int2str(year1(pos(i))) ' - ' int2str(year2(pos(i)))]); + vx = double(ncread(nc,'VX',[id1x id1y pos(i)],[id2x-id1x+1 id2y-id1y+1 1],[1 1 1])); + vy = double(ncread(nc,'VY',[id1x id1y pos(i)],[id2x-id1x+1 id2y-id1y+1 1],[1 1 1])); + vxdata(:,:,i) = permute(vx,[2 1 3]); + vydata(:,:,i) = permute(vy,[2 1 3]); +end +xdata=xdata(id1x:id2x); +ydata=ydata(id1y:id2y); + +disp([' -- Mouginot Time Series 1973 to 2018: interpolating']); +vxout = []; +vyout = []; +for i=1:length(pos), + disp([' step = ' int2str(i) '/' int2str(length(pos)) ', position = ' int2str(pos(i)) ', year = ' int2str(year1(pos(i))) ' - ' int2str(year2(pos(i)))]); + vxout = [vxout InterpFromGrid(xdata,ydata,vxdata(:,:,i),double(X),double(Y))]; + vyout = [vyout InterpFromGrid(xdata,ydata,vydata(:,:,i),double(X),double(Y))]; +end + +%return vel if only one output is requested +if nargout==1, + vxout = sqrt(vxout.^2+vyout.^2); +end Index: /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpPaolo2015.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpPaolo2015.m (revision 23873) +++ /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpPaolo2015.m (revision 23873) @@ -0,0 +1,187 @@ +function [dh_raw_out dh_fil_out T_out] = interpPaolo2015(X,Y,T,method) +%INTERPPAOLO2015 - interpolate observed (time series) height change [m] +% +% Inputs +% X,Y: spatial (scatter) coordinates +% T: time (see below the available years) +% ATTENTION: it is assumed that X and Y come in Polar Stereographic Projection (Std Latitude: 71S Meridian: 0E) +% +% Outputs +% dh_raw_out: interpolated raw height change at X,Y, for each time requested in T +% dh_fil_out: interpolated filtered height change at X,Y, for each time requested in T +% T_out: time related to dh_raw_out and dh_fil_out (see below) +% +% Available time series: +% +% 1 1994.038 +% 2 1994.285 +% 3 1994.534 +% 4 1994.786 +% 5 1995.038 +% 6 1995.285 +% 7 1995.534 +% 8 1995.786 +% 9 1996.038 +% 10 1996.287 +% 11 1996.536 +% 12 1996.787 +% 13 1997.038 +% 14 1997.285 +% 15 1997.534 +% 16 1997.786 +% 17 1998.038 +% 18 1998.285 +% 19 1998.534 +% 20 1998.786 +% 21 1999.038 +% 22 1999.285 +% 23 1999.534 +% 24 1999.786 +% 25 2000.038 +% 26 2000.287 +% 27 2000.536 +% 28 2000.787 +% 29 2001.038 +% 30 2001.285 +% 31 2001.534 +% 32 2001.786 +% 33 2002.038 +% 34 2002.285 +% 35 2002.534 +% 36 2002.786 +% 37 2003.038 +% 38 2003.285 +% 39 2003.534 +% 40 2003.786 +% 41 2004.038 +% 42 2004.287 +% 43 2004.536 +% 44 2004.787 +% 45 2005.038 +% 46 2005.285 +% 47 2005.534 +% 48 2005.786 +% 49 2006.038 +% 50 2006.285 +% 51 2006.534 +% 52 2006.786 +% 53 2007.038 +% 54 2007.285 +% 55 2007.534 +% 56 2007.786 +% 57 2008.038 +% 58 2008.287 +% 59 2008.536 +% 60 2008.787 +% 61 2009.038 +% 62 2009.285 +% 63 2009.534 +% 64 2009.786 +% 65 2010.038 +% 66 2010.285 +% 67 2010.534 +% 68 2010.786 +% 69 2011.038 +% 70 2011.285 +% 71 2011.534 +% 72 2011.786 +% +% +% Usage: +% % Get data at specific time: +% % In this example, T_out = [2006.038; 2007.038; 2008.038]. +% [dh_raw_out dh_fil_out T_out] = interpPaolo2015(md.mesh.x, md.mesh.y, [2006.038; 2007.038; 2008.038]); +% +% % Get all data in the provided years: +% % In this example, T_out = [2006.038; 2006.285; 2006.534; 2006.786; 2007.038; 2007.285; 2007.534; 2007.786]. +% [dh_raw dh_fil T_out] = interpPaolo2015(md.mesh.x, md.mesh.y, [2006; 2007]); +% +% % Get all data set: +% % In this example, T_out = [1994.038; ... ; 2011.786]. (all available time) +% [dh_raw dh_fil T_out] = interpPaolo2015(md.mesh.x, md.mesh.y); +% +% +% Info from ice_shelf_dh_v1.h5: +% The dataset is a rectangular grid (480 points in x, 80 points in y) with x- and y-axes being longitude and latitude, respectively. +% Longitude/latitude coordinates refer to the center of the grid cells. +% The grid has a resolution of lon x lat: 0.75 x 0.25 deg (~27 km at latitude -71). +% +% +% Data are (grids in HDF5, ice_shelf_dh_v1.h5): +% time : time coordinate [year; 72 values at 3-month time step] +% lon : x-coordinate [degrees east; range 0/360] +% lat : y-coordinate [degrees north; range -82/-62] +% height_raw : Raw time series of height change [m] +% height_filt : Filtered time series of height change [m] +% height_err : 2-standard-error time series [m] +% + +if nargin>4 | nargin<2, + error('nargin not supported yet!'); +end + +% read data +switch (oshostname()), + case {'ronne'} + h5 = '/home/ModelData/Antarctica/Paolo2015/ice_shelf_dh_v1.h5'; + otherwise + error('hostname not supported yet'); +end + +disp([' -- Paolo''s Time Series 1994 to 2012: loading data set']); +t_data = h5read(h5,'/time'); +lat_data = h5read(h5,'/lat'); +lon_data = h5read(h5,'/lon'); +dh_raw_data = h5read(h5,'/height_raw'); +dh_fil_data = h5read(h5,'/height_filt'); + +% set interpolation method +if nargin<4, + method = 'linear'; % default method +end + +% get the positions related to T +if nargin<3, + pos = 1:length(t_data); % all available data +else + % initial check %{{{ + if size(T,2)>1 | size(T,1)<1 | size(T,2)<1, + error('Size of input T not supported!'); + end %}}} + % Loop over T + pos = []; + epsilon = 5e-4; + for i=1:length(T), + % find specific time + flag = (T(i)-epsilont_data); + if ~any(flag), + % ok, find the time related to the requested year + flag = (T(i)==floor(t_data)); + end + if ~any(flag) + error(['requested time (' num2str(T(i)) ') not found in data set']) + end + pos = [pos; find(flag)]; + end + % Check if there is repeated positions + posunique = unique(pos); + if length(posunique)~=length(pos), + disp(' WARNING: found repeated positions in requested time'); + end +end + +% convert x/y to lat/lon: +[LAT, LON] = xy2ll(X,Y,-1); % attention: it is assumed that X and Y comes in Polar Stereographic Projection (Std Latitude: 71S Meridian: 0E) +posLON = find(LON<0); +LON(posLON) =360+LON(posLON); + +disp([' -- Paolo''s Time Series 1994 to 2012: interpolating in Lat/Long grid']); +dh_raw_out = []; +dh_fil_out = []; +for i=1:length(pos), + disp([' step = ' int2str(i) '/' int2str(length(pos)) ', position = ' int2str(pos(i)) ', year = ' num2str(t_data(pos(i)))]); + dh_raw_out = [dh_raw_out InterpFromGrid(lat_data(1,:),lon_data(:,1),dh_raw_data(:,:,pos(i)),LAT,LON,method)]; + dh_fil_out = [dh_fil_out InterpFromGrid(lat_data(1,:),lon_data(:,1),dh_fil_data(:,:,pos(i)),LAT,LON,method)]; +end + +T_out = t_data(pos); Index: /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpRACMO1km.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpRACMO1km.m (revision 23873) +++ /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpRACMO1km.m (revision 23873) @@ -0,0 +1,35 @@ +function [output] = interpRACMO1km(X,Y), + +switch oshostname(), + case {'ronne'} + rootname='/home/ModelData/Greenland/RACMO2_1km/SMB_MEAN1960-1989_150m.nc'; + otherwise + error('machine not supported yet'); +end +verbose = 1; + +xdata = double(ncread(rootname,'xaxis')); +ydata = double(ncread(rootname,'yaxis')); + +offset=2; + +xmin=min(X(:)); xmax=max(X(:)); +posx=find(xdata<=xmax); +id1x=max(1,find(xdata>=xmin,1)-offset); +id2x=min(numel(xdata),posx(end)+offset); + +ymin=min(Y(:)); ymax=max(Y(:)); +posy=find(ydata<=ymax); +id1y=max(1,find(ydata>=ymin,1)-offset); +id2y=min(numel(ydata),posy(end)+offset); + +if verbose, disp(' -- RACMO 1-km: reading smb'); end +data = double(ncread(rootname,'SMB',[id1x id1y],[id2x-id1x+1 id2y-id1y+1],[1 1]))'; +xdata=xdata(id1x:id2x); +ydata=ydata(id1y:id2y); +data(find(data==-9999))=NaN; + +if verbose, disp(' -- RACMO 1-km: interpolating (assuming rho_ice = 917 kg/m^3)'); end +%converting from mm / yr water eq to m/yr ice eq +data = data/1000 * 1000/917; +output = InterpFromGrid(xdata,ydata,data,double(X),double(Y)); Index: /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpRACMOant.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpRACMOant.m (revision 23873) +++ /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpRACMOant.m (revision 23873) @@ -0,0 +1,11 @@ +function smb = interpRACMOant(x,y); + + smbfile = '/home/ModelData/Antarctica/RACMO2SMB/SMB_RACMO2.3_1979_2011.nc'; + LAT=ncread(smbfile,'lat2d')'; + LON=ncread(smbfile,'lon2d')'; + SMB=ncread(smbfile,'SMB')'; + [X Y]=ll2xy(LAT,LON,-1,0,71); + + disp(' -- RACMO2.3 1979 - 2011: interpolating (assuming rho_ice = 917 kg/m^3)'); + rho_ice = 917; + smb = griddata(X,Y,SMB,x,y) / 917; Index: /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpREMA.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpREMA.m (revision 23873) +++ /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpREMA.m (revision 23873) @@ -0,0 +1,62 @@ +function sout = interpREMA(X,Y), + +switch oshostname(), + case {'ronne'} + remapath='/home/ModelData/Antarctica/REMA/REMA_200m_dem_filled.tif'; + otherwise + error('machine not supported yet'); +end + +usemap = 0; +if license('test','map_toolbox')==0, + disp('WARNING: map toolbox not installed, trying house code'); + usemap = 0; +elseif license('checkout','map_toolbox')==0 + disp('WARNING: map toolbox not available (checkout failed), trying house code'); + usemap = 0; +end + +if usemap, + [data,R] = geotiffread(remapath); + data=double(flipud(data)); + xdata=R.XLimWorld(1):R.DeltaX:R.XLimWorld(2); xdata=xdata(:); + xdata =(xdata(1:end-1)+xdata(2:end))/2; + ydata=R.YLimWorld(2):R.DeltaY:R.YLimWorld(1); ydata=flipud(ydata(:)); + ydata =(ydata(1:end-1)+ydata(2:end))/2; +else + + %Get image info + Tinfo = imfinfo(remapath); + N = Tinfo.Width; + M = Tinfo.Height; + dx = Tinfo.ModelPixelScaleTag(1); + dy = Tinfo.ModelPixelScaleTag(2); + minx = Tinfo.ModelTiepointTag(4); + maxy = Tinfo.ModelTiepointTag(5); + + %Generate vectors + xdata = minx + dx/2 + ((0:N-1).*dx); + ydata = maxy - dy/2 - ((M -1:-1:0).*dy); + ydata = fliplr(ydata); + + %Get pixels we are interested in + offset=2; + xmin=min(X(:)); xmax=max(X(:)); + posx=find(xdata<=xmax); + id1x=max(1,find(xdata>=xmin,1)-offset); + id2x=min(numel(xdata),posx(end)+offset); + + ymin=min(Y(:)); ymax=max(Y(:)); + posy=find(ydata>=ymin); + id1y=max(1,find(ydata<=ymax,1)-offset); + id2y=min(numel(ydata),posy(end)+offset); + + data = double(imread(remapath,'PixelRegion',{[id1y,id2y],[id1x,id2x]})); + xdata=xdata(id1x:id2x); + ydata=ydata(id1y:id2y); +end + +%convert no coverage data +data(find(data==-9999))=NaN; + +sout = InterpFromGrid(xdata,ydata,data,X,Y); Index: /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpRTopo2.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpRTopo2.m (revision 23873) +++ /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpRTopo2.m (revision 23873) @@ -0,0 +1,59 @@ +function [output] = interpRTopo2(X,Y,varargin), +%INTERPRTOPO2 - interp from RTOPO-2 onto X and Y +% +% Usage: +% bed = interpRTopo2(X,Y,varargin), +% +% varargin = 1 (Greenland), default +% -1 (Antarctica) + +switch oshostname(), + case {'ronne'} + rootname='/home/ModelData/Global/RTopo-2/RTopo-2.0.1_30sec_bedrock_topography.nc'; + otherwise + error('machine not supported yet'); +end +verbose = 1; + +if nargin==3, + hemisphere = varargin{1}; +else + hemisphere = +1; +end +if abs(hemisphere)~=1, + error('hemisphere should be +/-1'); +end + +if hemisphere==+1, + if verbose, disp(' -- RTopo-2: convert to lat/lon using Greenland projection'); end + [LAT, LON ] = xy2ll(double(X(:)),double(Y(:)),+1,45,70); +else + if verbose, disp(' -- RTopo-2: convert to lat/lon using Antarctica projection'); end + [LAT, LON ] = xy2ll(double(X(:)),double(Y(:)),-1,0,71); +end + +Y=reshape(LAT,size(X)); X=reshape(LON,size(X)); + +xdata = double(ncread(rootname,'lon')); +ydata = double(ncread(rootname,'lat')); + +offset=2; + +xmin=min(X(:)); xmax=max(X(:)); +posx=find(xdata<=xmax); +id1x=max(1,find(xdata>=xmin,1)-offset); +id2x=min(numel(xdata),posx(end)+offset); + +ymin=min(Y(:)); ymax=max(Y(:)); +posy=find(ydata<=ymax); +id1y=max(1,find(ydata>=ymin,1)-offset); +id2y=min(numel(ydata),posy(end)+offset); + +if verbose, disp(' -- RTopo-2: reading bed topography'); end +data = double(ncread(rootname,'bedrock_topography',[id1x id1y],[id2x-id1x+1 id2y-id1y+1],[1 1]))'; +xdata=xdata(id1x:id2x); +ydata=ydata(id1y:id2y); +data(find(data==-9999))=NaN; + +if verbose, disp(' -- RTopo-2: interpolating'); end +output = InterpFromGrid(xdata,ydata,data,double(X),double(Y)); Index: /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpRignot2012.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpRignot2012.m (revision 23873) +++ /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpRignot2012.m (revision 23873) @@ -0,0 +1,36 @@ +function [vxout vyout]= interpRignot2012(X,Y), + +filename = '/home/ModelData/Greenland/VelMouginot/RignotGreenland2012Vel.mat'; + + +%Figure out what subset of the matrix should be read +load(filename,'x','y'); +velfile = matfile(filename); + +offset=2; + +xmin=min(X(:)); xmax=max(X(:)); +posx=find(x<=xmax); +id1x=max(1,find(x>=xmin,1)-offset); +id2x=min(numel(x),posx(end)+offset); + +ymin=min(Y(:)); ymax=max(Y(:)); +%posy=find(y>=ymin); +%id1y=max(1,find(y<=ymax,1)-offset); +%id2y=min(numel(y),posy(end)+offset); +posy=find(y<=ymax); +id1y=max(1,find(y>=ymin,1)-offset); +id2y=min(numel(y),posy(end)+offset); + +vx = velfile.vx(id1y:id2y,id1x:id2x); +vy = velfile.vy(id1y:id2y,id1x:id2x); +x = x(id1x:id2x); +y = y(id1y:id2y); + +%load(filename); +vxout = InterpFromGrid(x,y,double(vx),X,Y); +vyout = InterpFromGrid(x,y,double(vy),X,Y); + +if nargout==1, + vxout = sqrt(vxout.^2+vyout.^2); +end Index: /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpRignotIceShelfMelt.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpRignotIceShelfMelt.m (revision 23873) +++ /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpRignotIceShelfMelt.m (revision 23873) @@ -0,0 +1,25 @@ +function output = interpRignotIceShelfMelt(X,Y,string), + +switch (oshostname()), + case {'ronne'} + rignotmelt='/home/ModelData/Antarctica/RignotMeltingrate/Ant_MeltingRate.nc'; + case {'thwaites','murdo','astrid'} + rignotmelt=['/home/seroussi/Data/Ant_MeltingRate.nc']; + otherwise + error('hostname not supported yet'); +end + +if nargin==2, + string = 'melt_actual'; +end + +disp([' -- Rignot Ice Shelf Melt: loading ' string]); +xdata = double(ncread(rignotmelt,'xaxis')); +ydata = double(ncread(rignotmelt,'yaxis')); + +disp([' -- Rignot Ice Shelf Melt: loading' string]); +data = double(ncread(rignotmelt,string))'; + +disp([' -- Rignot Ice Shelf Melt: interpolating ' string]); +output = InterpFromGrid(xdata,ydata,data,X(:),Y(:)); +output = reshape(output,size(X,1),size(X,2)); Index: /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpSeaRISE.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpSeaRISE.m (revision 23873) +++ /issm/trunk-jpl/src/m/contrib/morlighem/ronne/interpSeaRISE.m (revision 23873) @@ -0,0 +1,80 @@ +function [dataout] = interpSeaRISE(X,Y,string,varargin), +%INTERPSEARISE - interpolate SeaRISE data +% +% Available data: +% 1. sealeveltimes +% 2. dhdt +% 3. surfvelmag +% 4. balvelmag +% 5. oisotopestimes +% 6. bheatflx +% 7. presprcp +% 8. sealevel_time_series +% 9. usrf +% 10. mapping +% 11. surfvely +% 12. surfvelx +% 13. topg +% 14. landcover +% 15. temp_time_series +% 16. thk +% 17. time +% 18. oisotopes_time_series +% 19. runoff +% 20. smb +% 21. airtemp2m +% 22. surftemp +% +% Usage: +% [dataout] = interpBedmap2(X,Y,string) + +verbose=0; + +if nargin==3, + hemisphere = +1; +else + hemisphere = varargin{1}; +end + + +%read data +switch (oshostname()), + case {'ronne'} + if hemisphere==1, + searisenc='/home/ModelData/SeaRISE/Greenland_5km_dev1.2.nc'; + elseif hemisphere==-1, + searisenc='/home/ModelData/SeaRISE/Antarctica_5km_dev1.0.nc'; + end + case {'thwaites','murdo','astrid'} + if hemisphere==1, + searisenc='/u/astrid-r1b/ModelData/SeaRISE/Greenland5km_v1.2/Greenland_5km_dev1.2.nc'; + elseif hemisphere==-1, + searisenc='/u/astrid-r1b/ModelData/SeaRISE/Antarctica5km_shelves_v1.0/Antarctica_5km_dev1.0.nc'; + end + otherwise + error('hostname not supported yet'); +end + +%convert coordinates to SeaRISE projection +if verbose, disp(' -- SeaRISE: converting coordinates'); end +if hemisphere==1, + [LAT, LON ] = xy2ll(double(X(:)),double(Y(:)),+1,45,70); + [xproj,yproj] = ll2xy(LAT,LON ,+1,39,71); +elseif hemisphere==-1, + xproj=X; yproj=Y; +end + +if verbose, disp(' -- SeaRISE: loading coordinates'); end +xdata = double(ncread(searisenc,'x1'));%*1000; +ydata = double(ncread(searisenc,'y1'));%*1000; + +if verbose, disp([' -- SeaRISE: loading ' string]); end +data = double(ncread(searisenc,string))'; + +if verbose, disp([' -- SeaRISE: interpolating ' string]); end +if strcmpi(string,'LandMask'); + dataout = InterpFromGrid(xdata,ydata,data,xproj,yproj,'nearest'); +else + dataout = InterpFromGrid(xdata,ydata,data,xproj,yproj); +end +dataout = reshape(dataout,size(X,1),size(X,2));