Index: /issm/trunk-jpl/src/m/modeldata/InterpFromGrid.cpp =================================================================== --- /issm/trunk-jpl/src/m/modeldata/InterpFromGrid.cpp (revision 27548) +++ /issm/trunk-jpl/src/m/modeldata/InterpFromGrid.cpp (revision 27548) @@ -0,0 +1,567 @@ +/*Written by Mathieu Morlighem April 19th 2019*/ + +/*includes*/ +#include +#include +#include //for isnan +#include // for strcmp +#define f(m,n)\ + data[n*dataM+m] //Warning: matrix is transposed! + +/*Inputs{{{*/ +#define DATAX (mxArray*)prhs[0] +#define DATAY (mxArray*)prhs[1] +#define DATA (mxArray*)prhs[2] +#define INTERPX (mxArray*)prhs[3] +#define INTERPY (mxArray*)prhs[4] +#define METHOD (mxArray*)prhs[5] +/*}}}*/ +/*Outputs{{{*/ +#define INTERP (mxArray**)&plhs[0] +/*}}}*/ +/*threading structs{{{*/ +typedef struct{ + void* usr; + int my_thread; + int num_threads; +} pthread_handle; + +typedef struct{ + int dataM; + int dataN; + double* datax; + double* datay; + double* data; + int interpN; + double* interpx; + double* interpy; + double* interp; + int method; +} AppStruct; /*}}}*/ +/*Prototypes{{{*/ +void FetchMatrixPointer(double** pmatrix,int *pM,int *pN,const mxArray* dataref); +void FetchVectorPointer(double** pvector,int *pN,const mxArray* dataref); +void FetchString(char** pstring,const mxArray* dataref); +void WriteMatrix(mxArray** pdataref,double* matrix,int M,int N); +void WriteVector(mxArray** pdataref,double* vector,int N); +void* InterpFromGridt(void* vpthread_handle); +void LaunchThread(void* function(void*), void* usr,int num_threads); +bool binary_search_increasing(int* pindex,double target,double* list,int n); +bool binary_search_decreasing(int* pindex,double target,double* list,int n); +void dataderivatives(double* A,double* x,double* y,double* data,int M,int N, int m0, int m1,int m2,int m3, int n0, int n1,int n2,int n3); +/*}}}*/ + +void mexFunction(int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]){/*{{{*/ + + double *datax = NULL; + double *datay = NULL; + double *data = NULL; + int dataM,dataN; + double *interpx = NULL; + double *interpy = NULL; + double *interp = NULL; + int interpM,interpN; + + int num_threads = 20; + int test1,test2,test3,test4; + int method = 1; // 0 = nearest, 1 = bilinear, 2 = bicubic + + /*Check arguments to avoid crash*/ + if(nlhs>1 || (nrhs<5 || nrhs>6)) mexErrMsgTxt("Wrong usage"); + + /*Get variables from matlab to C*/ + FetchVectorPointer(&datax,&dataN,DATAX); + FetchVectorPointer(&datay,&dataM,DATAY); + FetchMatrixPointer(&data ,&test1,&test2,DATA); + FetchMatrixPointer(&interpx,&interpM,&interpN,INTERPX); + FetchMatrixPointer(&interpy,&test3,&test4,INTERPY); + if(!dataM*dataN) mexErrMsgTxt("data is empty"); + if(!interpM*interpN) mexErrMsgTxt("no interpolation requested"); + if(test1!=dataM) mexErrMsgTxt("x should have as many elements as there are columns in the data"); + if(test2!=dataN) mexErrMsgTxt("y should have as many elements as there are lines in the data"); + if(test3!=interpM) mexErrMsgTxt("interpolation locations (x,y) do not have the same size"); + if(test4!=interpN) mexErrMsgTxt("interpolation locations (x,y) do not have the same size"); + if(nrhs==6){ + char* method_string = NULL; + FetchString(&method_string,METHOD); + if(strcmp(method_string,"nearest")==0) method = 0; + else if(strcmp(method_string,"linear")==0) method = 1; + else if(strcmp(method_string,"cubic")==0) method = 2; + else{ + mexErrMsgTxt("Method not supported yet"); + } + + mxFree(method_string); + } + + /*Check inputs*/ + if(true){ + for(int i=0;imy_thread; + int num_threads = handle->num_threads; + + /*Recover struct*/ + AppStruct *usr = (AppStruct*)handle->usr; + int dataM = usr->dataM; + int dataN = usr->dataN; + double *datax = usr->datax; + double *datay = usr->datay; + double *data = usr->data; + int interpN = usr->interpN; + double *interpx = usr->interpx; + double *interpy = usr->interpy; + double *interp = usr->interp; + int method = usr->method; + + /*Intermediary*/ + double xprime,yprime; + double x,y,x0,x1,x2,x3,y0,y1,y2,y3; + double Q11,Q12; + double Q21,Q22; + double A[16]; + int m,n,m0,m1,m2,m3,n0,n1,n2,n3; + int oldm=-1,oldn=-1; + + /*Is our matrix inverted?*/ + bool invertx = (datax[1]-datax[0])<0 ? true:false; + bool inverty = (datay[1]-datay[0])<0 ? true:false; + + for(int idx=my_thread;idx=0 && n=0 && m (x1+x2)/2.){ + if(y > (y1+y2)/2.) + interp[idx] = f(m2,n2); + else + interp[idx] = f(m1,n2); + } + else{ + if(y > (y1+y2)/2.) + interp[idx] = f(m2,n1); + else + interp[idx] = f(m1,n1); + } + continue; + } + else if(method==1){ + /*Bilinear interpolation*/ + if(Q11==-9999 || Q12==-9999 || Q21==-9999 || Q22==-9999){ + interp[idx] = -9999; + continue; + } + + interp[idx] = + +f(m1,n1)*(x2-x)*(y2-y)/((x2-x1)*(y2-y1)) + +f(m1,n2)*(x-x1)*(y2-y)/((x2-x1)*(y2-y1)) + +f(m2,n1)*(x2-x)*(y-y1)/((x2-x1)*(y2-y1)) + +f(m2,n2)*(x-x1)*(y-y1)/((x2-x1)*(y2-y1)); + } + else{ + /*Bicubic interpolation*/ + if(invertx){n0=n+2; n3=n-1;} + else{ n0=n-1; n3=n+2; } + if(inverty){ m0=m+2; m3=m-1; } + else{ m0=m-1; m3=m+2; } + + if(n0<0 || n3>=dataN || m0<0 || m3>=dataM){ + interp[idx] = -9999.; + continue; + } + + /*Local coordinates (between 0 and 1)*/ + xprime = (x - datax[n1])/(datax[n2]-datax[n1]); + yprime = (y - datay[m1])/(datay[m2]-datay[m1]); + + /*Get derivatives at current pixel*/ + if(oldm!=m || oldn!=n){ + dataderivatives(&A[0],datax,datay,data,dataM,dataN,m0,m1,m2,m3,n0,n1,n2,n3); + oldm = m; + oldn = n; + } + + double a00 = A[0]; + double a10 = A[4]; + double a20 = -3*A[0]+3*A[1]-2*A[4]-A[5]; + double a30 = 2*A[0]-2*A[1]+A[4]+A[5]; + double a01 = A[8]; + double a11 = A[12]; + double a21 = -3*A[8]+3*A[9]-2*A[12]-A[13]; + double a31 = 2*A[8]-2*A[9]+A[12]+A[13]; + double a02 = -3*A[0]+3*A[2]-2*A[8]-A[10]; + double a12 = -3*A[4]+3*A[6]-2*A[12]-A[14]; + double a22 = 9*A[0]-9*A[1]-9*A[2]+9*A[3]+6*A[4]+3*A[5]-6*A[6]-3*A[7]+6*A[8]-6*A[9]+3*A[10]-3*A[11]+4*A[12]+2*A[13]+2*A[14]+A[15]; + double a32 =-6*A[0]+6*A[1]+6*A[2]-6*A[3]-3*A[4]-3*A[5]+3*A[6]+3*A[7]-4*A[8]+4*A[9]-2*A[10]+2*A[11]-2*A[12]-2*A[13]-A[14]-A[15]; + double a03 = 2*A[0]-2*A[2]+A[8]+A[10]; + double a13 = 2*A[4]-2*A[6]+A[12]+A[14]; + double a23 =-6*A[0]+6*A[1]+6*A[2]-6*A[3]-4*A[4]-2*A[5]+4*A[6]+2*A[7]-3*A[8]+3*A[9]-3*A[10]+3*A[11]-2*A[12]-A[13]-2*A[14]-A[15] ; + double a33 = 4*A[0]-4*A[1]-4*A[2]+4*A[3]+2*A[4]+2*A[5]-2*A[6]-2*A[7]+2*A[8]-2*A[9]+2*A[10]-2*A[11]+A[12]+A[13]+A[14]+A[15]; + + x1= xprime; + x2= x1*x1; + x3= x2*x1; + y1= yprime; + y2= y1*y1; + y3= y2*y1; + interp[idx] = (a00+a01*y1+a02*y2+a03*y3)+(a10+a11*y1+a12*y2+a13*y3)*x1+(a20+a21*y1+a22*y2+a23*y3)*x2+(a30+a31*y1+a32*y2+a33*y3)*x3; + } + } + else{ + interp[idx] = -9999.; + } + } + //if(my_thread==0) printf("\r interpolation progress = %5.1f%%\n",100.); + + return NULL; +}/*}}}*/ +/*binary_search_increasing {{{*/ +bool binary_search_increasing(int* pindex,double target,double* list,int n){ + + /*output*/ + int index; //index, if found + bool found=false; //found=0 if target is not found, 1 otherwise. + + /*intermediary*/ + int n0 = 0; + int n1 = int(n/2); + //int n1 = int((target-list[0])/(list[1]-list[0])); + int n2 = n-1; + + if(targetlist[n2]){ + found = true; + index = n; + } + else{ + while(!found){ + /*did we find the target?*/ + if(list[n1]<=target && list[n1+1]>=target){ + found = true; + index = n1; + break; + } + if(target < list[n1]){ + n2 = n1; + n1 = n0 + int((n2-n0)/2); + } + else{ + n0 = n1; + n1 = n0 + int((n2-n0)/2); + } + } + } + + /*Assign output pointers:*/ + *pindex=index; + + /*Return result: */ + return found; +}/*}}}*/ +/*binary_search_decreasing{{{*/ +bool binary_search_decreasing(int* pindex,double target,double* list,int n){ + + /*output*/ + int index; //index, if found + bool found=false; //found=0 if target is not found, 1 otherwise. + + /*intermediary*/ + int n0 = 0; + int n1 = int(n/2); + //int n1 = int((target-list[0])/(list[0]-list[1])); + int n2 = n-1; + + if (target>list[n0]){ + found = true; + index = -1; + } + else if(target=target && list[n1+1]<=target){ + found = true; + index = n1; + break; + } + if(target > list[n1]){ + n2 = n1; + n1 = n0 + int((n2-n0)/2); + } + else{ + n0 = n1; + n1 = n0 + int((n2-n0)/2); + } + } + } + + /*Assign output pointers:*/ + *pindex=index; + + /*Return result: */ + return found; +}/*}}}*/ +/*dataderivatives{{{*/ +void dataderivatives(double* A,double* x,double* y,double* data,int dataM,int dataN, + int m0, int m1,int m2,int m3, int n0, int n1,int n2,int n3){ + + /* i+1 + +-------+ f(1,1) + * | | | + * | |f(0,0) | + * i + +-------+ f(1,0) + * +--+-------+-----> x + * j j+1 + */ + + + /*Function at corners*/ + A[0] = f(m1,n1); // f(0,0) + A[1] = f(m1,n2); // f(1,0) + A[2] = f(m2,n1); // f(0,1) + A[3] = f(m2,n2); // f(1,1) + + /*x component of the gradient*/ + A[4] = .5*(f(m1,n2) - f(m1,n0));///(x[n2]-x[n0]); // dfdx(0,0) + A[5] = .5*(f(m1,n3) - f(m1,n1));///(x[n3]-x[n1]); // dfdx(1,0) + A[6] = .5*(f(m2,n2) - f(m2,n0));///(x[n2]-x[n0]); // dfdx(0,1) + A[7] = .5*(f(m2,n3) - f(m2,n1));///(x[n3]-x[n1]); // dfdx(1,1) + + /*y component of the gradient*/ + A[ 8] = .5*(f(m2,n1) - f(m0,n1));///(y[m2]-y[m0]); // dfdy(0,0) + A[ 9] = .5*(f(m2,n2) - f(m0,n2));///(y[m2]-y[m0]); // dfdy(1,0) + A[10] = .5*(f(m3,n1) - f(m1,n1));///(y[m3]-y[m1]); // dfdy(0,1) + A[11] = .5*(f(m3,n2) - f(m1,n2));///(y[m3]-y[m1]); // dfdy(1,1) + + /*cross-component of the gradient*/ + A[12] = .25*( (f(m2,n2) - f(m2,n0)) - (f(m0,n2) - f(m0,n0)) );///( (x[n2]-x[n0])*(y[m2]-y[m0]) ); // d2f/dxdy (0,0) + A[13] = .25*( (f(m2,n3) - f(m2,n1)) - (f(m0,n3) - f(m0,n1)) );///( (x[n3]-x[n1])*(y[m2]-y[m0]) ); // d2f/dxdy (1,0) + A[14] = .25*( (f(m3,n2) - f(m3,n0)) - (f(m1,n2) - f(m1,n0)) );///( (x[n2]-x[n0])*(y[m3]-y[m1]) ); // d2f/dxdy (0,1) + A[15] = .25*( (f(m3,n3) - f(m3,n1)) - (f(m1,n3) - f(m1,n1)) );///( (x[n3]-x[n1])*(y[m3]-y[m1]) ); // d2f/dxdy (1,1) +}/*}}}*/ +/*LaunchThread{{{*/ +void LaunchThread(void* function(void*), void* usr,int num_threads){ + + int i; + int *status = NULL; + pthread_t *threads = NULL; + pthread_handle *handles = NULL; + + /*dynamically allocate: */ + threads=(pthread_t*)mxMalloc(num_threads*sizeof(pthread_t)); + handles=(pthread_handle*)mxMalloc(num_threads*sizeof(pthread_handle)); + + for(i=0;i=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(path,'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/modeldata/interpBamber2001.m =================================================================== --- /issm/trunk-jpl/src/m/modeldata/interpBamber2001.m (revision 27548) +++ /issm/trunk-jpl/src/m/modeldata/interpBamber2001.m (revision 27548) @@ -0,0 +1,36 @@ +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'; + case {'totten'} + bamber2001bedpath ='/totten_1/ModelData/Greenland/Bamber2001/bedrock.mat'; + bamber2001thxpath ='/totten_1/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/modeldata/interpBamber2009.m =================================================================== --- /issm/trunk-jpl/src/m/modeldata/interpBamber2009.m (revision 27548) +++ /issm/trunk-jpl/src/m/modeldata/interpBamber2009.m (revision 27548) @@ -0,0 +1,21 @@ +function demout = interpBamber2009(X, Y) +%INTERPBAMBER2009 - interpolate surface dem of Bamber 2009 +% +% Surface dem Nominal year 2004 (WGS84, no firn correction) +% +% Usage: +% demout = interpBamber2009(X, Y) + +switch oshostname(), + case {'totten'} + bamber2009path ='/totten_1/ModelData/Antarctica/Bamber2009DEM/krigged_dem_nsidc.mat'; + otherwise + error('machine not supported yet'); +end + +%Convert to Bamber's projections +%disp(' -- Bamber2009: loading dem'); +load(bamber2009path); + +disp(' -- Bamber2009: interpolating dem (WGS84)'); +demout = InterpFromGrid(x, y, surfacedem, X, Y); Index: /issm/trunk-jpl/src/m/modeldata/interpBamber2013.m =================================================================== --- /issm/trunk-jpl/src/m/modeldata/interpBamber2013.m (revision 27548) +++ /issm/trunk-jpl/src/m/modeldata/interpBamber2013.m (revision 27548) @@ -0,0 +1,52 @@ +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'; + case {'totten'} + bamber2013nc='/totten_1/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/modeldata/interpBedmachineAntarctica.m =================================================================== --- /issm/trunk-jpl/src/m/modeldata/interpBedmachineAntarctica.m (revision 27548) +++ /issm/trunk-jpl/src/m/modeldata/interpBedmachineAntarctica.m (revision 27548) @@ -0,0 +1,165 @@ +function output = interpBedmachineAntarctica(X,Y,string,method,ncdate) +%INTERPBEDMACHINEANTARCTICA - interpolate BedMachine data onto X and Y +% +% Examples: +% bed = interpBedmachineAntarctica(X,Y,'bed'); +% surface = interpBedmachineAntarctica(X,Y,'surface'); +% thickness = interpBedmachineAntarctica(X,Y,'thickness'); +% mask = interpBedmachineAntarctica(X,Y,'mask'); +% mask = interpBedmachineAntarctica(X,Y,'mask','nearest','../Data/BedMachineAntarctica_2020-07-15_v02.nc'); +% +% - mask: 0 ocean, 1 land (ice free), 2 grounded ice, 3 floating ice +% - source: 1 IBCSO/RTopo-2, 2 MC, 3 interpolation, 4 hydrostatic eq, +% 5 Streamline diffusion, 6 Gravity inversion +% - optional 4th input argument: interpolation method. +% Supported interpolation methos: 'linear','cubic','nearest' +% - optional 5th input argument: path to dataset. +% +% Version 11/30/2018 Mathieu Morlighem mmorligh@uci.edu + +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='2020-07-15'; %BedMachine v2 + ncdate='v3.5'; %Official v3 release +end +basename = 'BedMachineAntarctica'; + +if nargin==5 + ncfile = ncdate; +else + %List of common paths to try + paths = {... + ['/u/astrid-r1b/ModelData/BedMachine/' basename '-' ncdate '.nc'],... + ['/home/ModelData/Antarctica/BedMachine/' basename '-' ncdate '.nc'],... + ['/totten_1/ModelData/Antarctica/BedMachine/' basename '-' ncdate '.nc'],... + ['/Users/larour/ModelData/BedMachine/' basename '-' ncdate '.nc'],... + ['./' basename '-' ncdate '.nc'],... + }; + + found = 0; + for i=1:numel(paths) + if exist(paths{i},'file') + ncfile = paths{i}; + found = 1; + break; + end + end + + if ~found + error(['Could not find ' basename '-' ncdate '.nc, you can add the path to the list or provide its path as a 5th argument']); + end +end + +disp([' -- BedMachine Antarctica version: ' ncdate]); +xdata = double(ncread(ncfile,'x')); +ydata = double(ncread(ncfile,'y')); + +offset=2; + +xmin=min(X(:)); xmax=max(X(:)); +posx=find(xdata<=xmax); +if isempty(posx), posx=numel(xdata); end +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); +if isempty(posy), posy=numel(ydata); end +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(ncfile,'mask'))'; + data = double(ncread(ncfile,'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(ncfile,string))'; + data = double(ncread(ncfile,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 + output = InterpFromGrid(xdata,ydata,data,double(X),double(Y),'nearest'); + %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(ncfile,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]); +disp([' -- Interpolation method: ' method]); +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 + +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); + 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; + +%For Eric's computer (using Binary files) +elseif exist('/Users/larour/ModelData/BedMap2/bedmap2_bin/','dir') + % ================================ OLD =============================================== + path='/Users/larour/ModelData/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 +else + error('not supported'); +end Index: /issm/trunk-jpl/src/m/modeldata/interpChuter2015.m =================================================================== --- /issm/trunk-jpl/src/m/modeldata/interpChuter2015.m (revision 27548) +++ /issm/trunk-jpl/src/m/modeldata/interpChuter2015.m (revision 27548) @@ -0,0 +1,13 @@ +function output = interpChuter2015(X,Y), + +ncfile='/home/ModelData/Antarctica/ChuterBamberIceShelfH/ChuterBamber_2015_CS2_ice_equivalent_ice_shelf_thickness_Rignot_gl.nc'; +verbose = 0; + +if verbose, disp(' -- Chuter2015: loading coordinates'); end +xdata = double(ncread(ncfile,'x_dimensions'))'; +ydata = double(ncread(ncfile,'y_dimensions'))'; + +if verbose, disp([' -- Chuter2015: loading thickenss']); end +data = double(ncread(ncfile,'ice_shelf_thickness'))'; +if verbose, disp([' -- Chuter2015: interpolating ' string]); end +output = InterpFromGrid(xdata(1,:),ydata(:,1),data,X,Y); Index: /issm/trunk-jpl/src/m/modeldata/interpDTU19MDT.m =================================================================== --- /issm/trunk-jpl/src/m/modeldata/interpDTU19MDT.m (revision 27548) +++ /issm/trunk-jpl/src/m/modeldata/interpDTU19MDT.m (revision 27548) @@ -0,0 +1,35 @@ +function mdt = interpDTU19MDT(X,Y,varargin); + +switch oshostname(), + case {'ronne'} + rootname='/ronne_2/home/ModelData/Global/DTU19MDT/dtu19mdt.mat'; + case {'totten'} + rootname='/totten_1/ModelData/Global/DTU19MDT/dtu19mdt.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(' -- DTU19MDT: convert to lat/lon using Greenland projection'); end + [LAT, LON ] = xy2ll(double(X(:)),double(Y(:)),+1,45,70); +else + if verbose, disp(' -- DTU19MDT: 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); +LAT=reshape(LAT,size(X)); +LON=reshape(LON,size(X)); + +if verbose, disp(' -- DTU19MDT: loading DTU19MDT'); end +A=load(rootname); + +if verbose, disp(' -- DTU19MDT: interpolating'); end +mdt = InterpFromGrid(A.lon_ext,A.lat_ext,A.mdt_ext,LON,LAT); Index: /issm/trunk-jpl/src/m/modeldata/interpDhdt.m =================================================================== --- /issm/trunk-jpl/src/m/modeldata/interpDhdt.m (revision 27548) +++ /issm/trunk-jpl/src/m/modeldata/interpDhdt.m (revision 27548) @@ -0,0 +1,55 @@ +function out = interpDhdt(X,Y), + +switch oshostname(), + case {'ronne'} + dhdtpath='/home/ModelData/Greenland/DHDT/dhdt0306.tif'; + case {'totten'} + dhdtpath='/totten_1/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/modeldata/interpFromGeotiff.m =================================================================== --- /issm/trunk-jpl/src/m/modeldata/interpFromGeotiff.m (revision 27548) +++ /issm/trunk-jpl/src/m/modeldata/interpFromGeotiff.m (revision 27548) @@ -0,0 +1,70 @@ +function dataout = interpFromGeotiff(geotiffname,X,Y,nanValue), + +if nargin < 4 + nanValue = 10^30; +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(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(1).Width; + M = Tinfo(1).Height; + dx = Tinfo(1).ModelPixelScaleTag(1); + dy = Tinfo(1).ModelPixelScaleTag(2); + minx = Tinfo(1).ModelTiepointTag(4); + maxy = Tinfo(1).ModelTiepointTag(5); + + %Generate vectors + xdata = minx + dx/2 + ((0:N-1).*dx); + ydata = maxy - dy/2 - ((M -1:-1:0).*dy); + + %Read image + if 1 + assert(dx>0); assert(dy>0); + 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(geotiffname,'PixelRegion',{[id1y,id2y],[id1x,id2x]})); + xdata=xdata(id1x:id2x); + ydata=ydata(id1y:id2y); + else + data=double(flipud(imread(geotiffname))); + end + if nanValue > 0 + data(find(abs(data)>=nanValue))=NaN; + else + data(find(data<=nanValue))=NaN; + end + +end + +dataout = InterpFromGrid(xdata,ydata,data,X,Y); +dataout(dataout==-9999)=NaN; Index: /issm/trunk-jpl/src/m/modeldata/interpGeoid.m =================================================================== --- /issm/trunk-jpl/src/m/modeldata/interpGeoid.m (revision 27548) +++ /issm/trunk-jpl/src/m/modeldata/interpGeoid.m (revision 27548) @@ -0,0 +1,36 @@ +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'; + case {'totten'} + rootname='/totten_1/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/modeldata/interpGimpdem.m =================================================================== --- /issm/trunk-jpl/src/m/modeldata/interpGimpdem.m (revision 27548) +++ /issm/trunk-jpl/src/m/modeldata/interpGimpdem.m (revision 27548) @@ -0,0 +1,71 @@ +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'; + case {'totten'} + %howatpath='/totten_1/ModelData/Greenland/gimpdem/gimpdem_90m_v01.1.tif'; + howatpath='/totten_1/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/modeldata/interpGimpicemask.m =================================================================== --- /issm/trunk-jpl/src/m/modeldata/interpGimpicemask.m (revision 27548) +++ /issm/trunk-jpl/src/m/modeldata/interpGimpicemask.m (revision 27548) @@ -0,0 +1,71 @@ +function sout = interpGimpicemask(X,Y), + +switch oshostname(), + case {'ronne'} + howatpath='/home/ModelData/Greenland/gimpmask/GimpIceMask_90m.tif'; + case {'totten'} + howatpath='/totten_1/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/modeldata/interpGimpoceanmask.m =================================================================== --- /issm/trunk-jpl/src/m/modeldata/interpGimpoceanmask.m (revision 27548) +++ /issm/trunk-jpl/src/m/modeldata/interpGimpoceanmask.m (revision 27548) @@ -0,0 +1,71 @@ +function sout = interpGimpoceanmask(X,Y), + +switch oshostname(), + case {'ronne'} + howatpath='/home/ModelData/Greenland/gimpmask/GimpOceanMask_90m.tif'; + case {'totten'} + howatpath='/totten_1/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/modeldata/interpGrIMP.m =================================================================== --- /issm/trunk-jpl/src/m/modeldata/interpGrIMP.m (revision 27548) +++ /issm/trunk-jpl/src/m/modeldata/interpGrIMP.m (revision 27548) @@ -0,0 +1,66 @@ +function sout = interpGimpdem(X,Y), + +switch oshostname(), + case {'totten'} + howatpath='/totten_1/ModelData/Greenland/GrIMP/GrIMP_100m.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/modeldata/interpGridsCReSIS.m =================================================================== --- /issm/trunk-jpl/src/m/modeldata/interpGridsCReSIS.m (revision 27548) +++ /issm/trunk-jpl/src/m/modeldata/interpGridsCReSIS.m (revision 27548) @@ -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/modeldata/interpGriggs2013.m =================================================================== --- /issm/trunk-jpl/src/m/modeldata/interpGriggs2013.m (revision 27548) +++ /issm/trunk-jpl/src/m/modeldata/interpGriggs2013.m (revision 27548) @@ -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/modeldata/interpIBCSO.m =================================================================== --- /issm/trunk-jpl/src/m/modeldata/interpIBCSO.m (revision 27548) +++ /issm/trunk-jpl/src/m/modeldata/interpIBCSO.m (revision 27548) @@ -0,0 +1,41 @@ +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'; + case {'totten'} + ncpath='/totten_1/ModelData/Antarctica/IBCSO/ibcso_v1_bed.grd'; + sidpath='/totten_1/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/modeldata/interpIBCSO2.m =================================================================== --- /issm/trunk-jpl/src/m/modeldata/interpIBCSO2.m (revision 27548) +++ /issm/trunk-jpl/src/m/modeldata/interpIBCSO2.m (revision 27548) @@ -0,0 +1,27 @@ +function [bedout sid] = interpIBCSO2(X,Y), + +%read data +switch (oshostname()), + case {'totten'} + ncpath='/totten_1/ModelData/Antarctica/IBCSO2/IBCSO_v2_bed.nc'; + sidpath='/totten_1/ModelData/Antarctica/IBCSO2/IBCSO_v2_TID.nc'; + otherwise + error('hostname not supported yet'); +end + +disp(' -- IBCSOv2: Changing Coordinate system from 3031 to 9354'); +[X Y]=CoordTransform(double(X),double(Y),'EPSG:3031','EPSG:9354'); + +disp(' -- IBCSOv2: loading bathymetry'); +xdata = double(ncread(ncpath,'x')); +ydata = double(ncread(ncpath,'y')); +data = double(ncread(ncpath,'z'))'; +disp(' -- IBCSOv2: interpolating bed'); +bedout = InterpFromGrid(xdata,ydata,data,double(X),double(Y)); + +if nargout==2, + disp(' -- IBCSOv2: bathymetry sid'); + data = ncread(sidpath,'tid')'; + disp(' -- IBCSOv2: interpolating sids'); + sid = InterpFromGrid(xdata,ydata,data,double(X),double(Y),'nearest'); +end Index: /issm/trunk-jpl/src/m/modeldata/interpJakobsson2012.m =================================================================== --- /issm/trunk-jpl/src/m/modeldata/interpJakobsson2012.m (revision 27548) +++ /issm/trunk-jpl/src/m/modeldata/interpJakobsson2012.m (revision 27548) @@ -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/modeldata/interpJakobsson2020.m =================================================================== --- /issm/trunk-jpl/src/m/modeldata/interpJakobsson2020.m (revision 27548) +++ /issm/trunk-jpl/src/m/modeldata/interpJakobsson2020.m (revision 27548) @@ -0,0 +1,49 @@ +function [bedout sourceout] = interpJakobsson2020(X,Y,string), + +switch oshostname(), + case {'ronne'} + ncpath ='/home/ModelData/Greenland/IBCAO/IBCAO_v4_200m.nc'; + case {'totten'} + ncpath ='/totten_1/ModelData/Greenland/IBCAO/IBCAO_v4_200m.nc'; + otherwise + error('machine not supported yet'); +end + +%Convert to IBCAO projections +disp(' -- Jakobsson2020: converting coordinates'); +[LAT, LON ] = xy2ll(double(X(:)),double(Y(:)),+1,45,70); +[x0075,y0075] = ll2xy(LAT,LON,+1,0,75); + +disp(' -- Jakobsson2020: loading coordinates'); +xdata = double(ncread(ncpath,'x')); +ydata = double(ncread(ncpath,'y')); + +offset=2; + +xmin=min(x0075(:)); xmax=max(x0075(:)); +posx=find(xdata<=xmax); +id1x=max(1,find(xdata>=xmin,1)-offset); +id2x=min(numel(xdata),posx(end)+offset); + +ymin=min(y0075(:)); ymax=max(y0075(:)); +posy=find(ydata>=ymin); +id1y=max(1,find(ydata<=ymax,1)-offset); +id2y=min(numel(ydata),posy(end)+offset); + +disp([' -- Jakobsson2020: loading bathymetry']); +data = double(ncread(ncpath,'z',[id1x id1y],[id2x-id1x+1 id2y-id1y+1],[1 1]))'; +xdata=xdata(id1x:id2x); +ydata=ydata(id1y:id2y); + +disp(' -- Jakobsson2020: interpolating bed'); +bedout = InterpFromGrid(xdata,ydata,data,x0075,y0075); +bedout = reshape(bedout,size(X,1),size(X,2)); + +if nargout==2, + ncpath ='/totten_1/ModelData/Greenland/IBCAO/IBCAO_v4_200m_TID.nc'; + disp(' -- Jakobsson2020: loading source'); + data = double(ncread(ncpath,'z',[id1x id1y],[id2x-id1x+1 id2y-id1y+1],[1 1]))'; + disp(' -- Jakobsson2020: 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/modeldata/interpJoughin.m =================================================================== --- /issm/trunk-jpl/src/m/modeldata/interpJoughin.m (revision 27548) +++ /issm/trunk-jpl/src/m/modeldata/interpJoughin.m (revision 27548) @@ -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/modeldata/interpJoughinCompositeGreenland.m =================================================================== --- /issm/trunk-jpl/src/m/modeldata/interpJoughinCompositeGreenland.m (revision 27548) +++ /issm/trunk-jpl/src/m/modeldata/interpJoughinCompositeGreenland.m (revision 27548) @@ -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 = '/totten_1/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/modeldata/interpJoughinMosaic.m =================================================================== --- /issm/trunk-jpl/src/m/modeldata/interpJoughinMosaic.m (revision 27548) +++ /issm/trunk-jpl/src/m/modeldata/interpJoughinMosaic.m (revision 27548) @@ -0,0 +1,36 @@ +function [vxout vyout] = interpJoughinMosaic(X,Y), + +switch oshostname(), + case {'ronne'} + filename = '/home/ModelData/Greenland/VelJoughin/IanGreenVel.mat'; + case {'totten'} + filename = '/totten_1/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/modeldata/interpMartos2017.m =================================================================== --- /issm/trunk-jpl/src/m/modeldata/interpMartos2017.m (revision 27548) +++ /issm/trunk-jpl/src/m/modeldata/interpMartos2017.m (revision 27548) @@ -0,0 +1,18 @@ +function out = interpMartos2017(X,Y) +%INTERPMARTOS2017 - interpolate geothermal heat flux +% +% Usage: +% out = interpMartos2017(X,Y) + +switch oshostname(), + case {'ronne'} + gtfpath='/home/ModelData/Antarctica/GeothermalMartos/Antarctic_GHF.xyz'; + otherwise + error('machine not supported yet'); +end + +%Load data +data = load(gtfpath); + +%Interpolate using nearest neighbor (dataset stops at ocean boundary!) +out = Kriging(data(:,1),data(:,2),data(:,3),X,Y,'output','nearestneighbor')/1e3; %from mW/m2 to W/m2 Index: /issm/trunk-jpl/src/m/modeldata/interpMouginotAnt2016.m =================================================================== --- /issm/trunk-jpl/src/m/modeldata/interpMouginotAnt2016.m (revision 27548) +++ /issm/trunk-jpl/src/m/modeldata/interpMouginotAnt2016.m (revision 27548) @@ -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/modeldata/interpMouginotAnt2017.m =================================================================== --- /issm/trunk-jpl/src/m/modeldata/interpMouginotAnt2017.m (revision 27548) +++ /issm/trunk-jpl/src/m/modeldata/interpMouginotAnt2017.m (revision 27548) @@ -0,0 +1,41 @@ +function [vxout vyout]= interpMouginotAnt2017(X,Y), + +%read data +switch (oshostname()), + case {'ronne'} + nc = '/home/ModelData/Antarctica/MouginotVel/vel_nsidc.CF16_2.nc'; + case {'totten'} + nc = '/totten_1/ModelData/Antarctica/MouginotVel/vel_nsidc.CF16_2.nc'; + otherwise + error('hostname not supported yet'); +end + +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); + +disp([' -- Mouginot 2017: loading velocities']); +vxdata = double(ncread(nc,'VX',[id1x id1y],[id2x-id1x+1 id2y-id1y+1],[1 1]))'; +vydata = double(ncread(nc,'VY',[id1x id1y],[id2x-id1x+1 id2y-id1y+1],[1 1]))'; +xdata=xdata(id1x:id2x); +ydata=ydata(id1y:id2y); + +disp([' -- Mouginot 2017: interpolating ']); +vxout = InterpFromGrid(xdata,ydata,vxdata,double(X),double(Y)); +vyout = InterpFromGrid(xdata,ydata,vydata,double(X),double(Y)); + +%return vel if only one output is requested +if nargout==1, + vxout = sqrt(vxout.^2+vyout.^2); +end Index: /issm/trunk-jpl/src/m/modeldata/interpMouginotAnt2019.m =================================================================== --- /issm/trunk-jpl/src/m/modeldata/interpMouginotAnt2019.m (revision 27548) +++ /issm/trunk-jpl/src/m/modeldata/interpMouginotAnt2019.m (revision 27548) @@ -0,0 +1,41 @@ +function [vxout vyout]= interpMouginotAnt2019(X,Y), + +%read data +switch (oshostname()), + case {'ronne'} + nc = '/home/ModelData/Antarctica/MouginotVel/v_mix.v13Mar2019.nc'; + case {'totten'} + nc = '/totten_1/ModelData/Antarctica/MouginotVel/v_mix.v8Jul2019.nc'; + otherwise + error('hostname not supported yet'); +end + +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); + +disp([' -- Mouginot 2019: loading velocities']); +vxdata = double(ncread(nc,'VX',[id1x id1y],[id2x-id1x+1 id2y-id1y+1],[1 1]))'; +vydata = double(ncread(nc,'VY',[id1x id1y],[id2x-id1x+1 id2y-id1y+1],[1 1]))'; +xdata=xdata(id1x:id2x); +ydata=ydata(id1y:id2y); + +disp([' -- Mouginot 2019: interpolating ']); +vxout = InterpFromGrid(xdata,ydata,vxdata,double(X),double(Y)); +vyout = InterpFromGrid(xdata,ydata,vydata,double(X),double(Y)); + +%return vel if only one output is requested +if nargout==1, + vxout = sqrt(vxout.^2+vyout.^2); +end Index: /issm/trunk-jpl/src/m/modeldata/interpMouginotAntTimeSeries1973to2018.m =================================================================== --- /issm/trunk-jpl/src/m/modeldata/interpMouginotAntTimeSeries1973to2018.m (revision 27548) +++ /issm/trunk-jpl/src/m/modeldata/interpMouginotAntTimeSeries1973to2018.m (revision 27548) @@ -0,0 +1,196 @@ +function [vxout vyout errxout erryout stdxout stdyout]= 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 examples: +% [vxout vyout]= interpMouginotAntTimeSeries1973to2018(md.mesh.x,md.mesh.y,[1973 1975; 1973 1988; 1991 1992; 2011 2012]); +% [vel]= interpMouginotAntTimeSeries1973to2018(md.mesh.x,md.mesh.y,[1986; 1991; 1995; 2000]); +% [vxout vyout errxout erryout stdxout stdyout]= interpMouginotAntTimeSeries1973to2018(md.mesh.x,md.mesh.y,[1986; 1991; 1995; 2000]); + +%read data +switch (oshostname()), + case {'ronne'} + nc = '/home/ModelData/Antarctica/MouginotVel/ASE_TimeSeries_1973-2018.nc'; + case {'totten'} + nc = '/totten_1/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 +if nargout~=1 & nargout~=2 & nargout~=6 + error('nargout not supported!'); +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 = []; +if nargout==6 % it includes ERRX, ERRY, STDX and STDY + errxdata = []; + errydata = []; + stdxdata = []; + stdydata = []; +end +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]); + if nargout==6 % it includes ERRX, ERRY, STDX and STDY + errx = double(ncread(nc,'ERRX',[id1x id1y pos(i)],[id2x-id1x+1 id2y-id1y+1 1],[1 1 1])); + erry = double(ncread(nc,'ERRY',[id1x id1y pos(i)],[id2x-id1x+1 id2y-id1y+1 1],[1 1 1])); + stdx = double(ncread(nc,'STDX',[id1x id1y pos(i)],[id2x-id1x+1 id2y-id1y+1 1],[1 1 1])); + stdy = double(ncread(nc,'STDY',[id1x id1y pos(i)],[id2x-id1x+1 id2y-id1y+1 1],[1 1 1])); + errxdata(:,:,i) = permute(errx,[2 1 3]); + errydata(:,:,i) = permute(erry,[2 1 3]); + stdxdata(:,:,i) = permute(stdx,[2 1 3]); + stdydata(:,:,i) = permute(stdy,[2 1 3]); + end +end +xdata=xdata(id1x:id2x); +ydata=ydata(id1y:id2y); + +disp([' -- Mouginot Time Series 1973 to 2018: interpolating']); +vxout = []; +vyout = []; +if nargout==6 % it includes ERRX, ERRY, STDX and STDY + errxout = []; + erryout = []; + stdxout = []; + stdyout = []; +end +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))]; + if nargout==6 % it includes ERRX, ERRY, STDX and STDY + errxout = [errxout InterpFromGrid(xdata,ydata,errxdata(:,:,i),double(X),double(Y))]; + erryout = [erryout InterpFromGrid(xdata,ydata,errydata(:,:,i),double(X),double(Y))]; + stdxout = [stdxout InterpFromGrid(xdata,ydata,stdxdata(:,:,i),double(X),double(Y))]; + stdyout = [stdyout InterpFromGrid(xdata,ydata,stdydata(:,:,i),double(X),double(Y))]; + end +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/modeldata/interpPaolo2015.m =================================================================== --- /issm/trunk-jpl/src/m/modeldata/interpPaolo2015.m (revision 27548) +++ /issm/trunk-jpl/src/m/modeldata/interpPaolo2015.m (revision 27548) @@ -0,0 +1,193 @@ +function [dh_raw_out dh_fil_out T_out] = interpPaolo2015(X,Y,T,method) +%INTERPPAOLO2015 - interpolate observed (time series) height change [m] +% +% Time series are average height changes [m] with respect to 1994 every three months (72 time steps) +% +% 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 + if size(X,1)>1 & size(X,2)>1 + error('Size of input X not supported! X and Y should be vectors'); + 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/modeldata/interpRACMO1km.m =================================================================== --- /issm/trunk-jpl/src/m/modeldata/interpRACMO1km.m (revision 27548) +++ /issm/trunk-jpl/src/m/modeldata/interpRACMO1km.m (revision 27548) @@ -0,0 +1,37 @@ +function [output] = interpRACMO1km(X,Y), + +switch oshostname(), + case {'ronne'} + rootname='/home/ModelData/Greenland/RACMO2_1km/SMB_MEAN1960-1989_150m.nc'; + case {'totten'} + rootname='/totten_1/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/modeldata/interpRACMOant.m =================================================================== --- /issm/trunk-jpl/src/m/modeldata/interpRACMOant.m (revision 27548) +++ /issm/trunk-jpl/src/m/modeldata/interpRACMOant.m (revision 27548) @@ -0,0 +1,18 @@ +function smb = interpRACMOant(x,y); + +switch oshostname(), + case {'ronne'} + smbfile = '/home/ModelData/Antarctica/RACMO2SMB/SMB_RACMO2.3_1979_2011.nc'; + case {'totten'} + smbfile = '/totten_1/ModelData/Antarctica/RACMO2SMB/SMB_RACMO2.3_1979_2011.nc'; + otherwise + error('machine not supported yet'); +end + 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/modeldata/interpREMA.m =================================================================== --- /issm/trunk-jpl/src/m/modeldata/interpREMA.m (revision 27548) +++ /issm/trunk-jpl/src/m/modeldata/interpREMA.m (revision 27548) @@ -0,0 +1,64 @@ +function sout = interpREMA(X,Y), + +switch oshostname(), + case {'ronne'} + remapath='/home/ModelData/Antarctica/REMA/REMA_200m_dem_filled.tif'; + case {'totten'} + remapath='/totten_1/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/modeldata/interpRTopo2.m =================================================================== --- /issm/trunk-jpl/src/m/modeldata/interpRTopo2.m (revision 27548) +++ /issm/trunk-jpl/src/m/modeldata/interpRTopo2.m (revision 27548) @@ -0,0 +1,61 @@ +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'; + case {'totten'} + rootname='/totten_1/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/modeldata/interpRignot2012.m =================================================================== --- /issm/trunk-jpl/src/m/modeldata/interpRignot2012.m (revision 27548) +++ /issm/trunk-jpl/src/m/modeldata/interpRignot2012.m (revision 27548) @@ -0,0 +1,36 @@ +function [vxout vyout]= interpRignot2012(X,Y), + +filename = '/totten_1/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/modeldata/interpRignotIceShelfMelt.m =================================================================== --- /issm/trunk-jpl/src/m/modeldata/interpRignotIceShelfMelt.m (revision 27548) +++ /issm/trunk-jpl/src/m/modeldata/interpRignotIceShelfMelt.m (revision 27548) @@ -0,0 +1,31 @@ +function output = interpRignotIceShelfMelt(X,Y,string) +%INTERPRIGNOTICESHELFMELT - interp melt rates from Rignot et al. 2013 +% +% Usage: +% output = interpRignotIceShelfMelt(X,Y) + +switch (oshostname()) + case {'ronne'} + rignotmelt='/home/ModelData/Antarctica/RignotMeltingrate/Ant_MeltingRate.nc'; + case {'totten'} + rignotmelt='/totten_1/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/modeldata/interpSeaRISE.m =================================================================== --- /issm/trunk-jpl/src/m/modeldata/interpSeaRISE.m (revision 27548) +++ /issm/trunk-jpl/src/m/modeldata/interpSeaRISE.m (revision 27548) @@ -0,0 +1,86 @@ +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 + case {'totten'} + if hemisphere==1, + searisenc='/totten_1/ModelData/SeaRISE/Greenland_5km_dev1.2.nc'; + elseif hemisphere==-1, + searisenc='/totten_1/ModelData/SeaRISE/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)); Index: /issm/trunk-jpl/src/m/modeldata/interpShepherd2019.m =================================================================== --- /issm/trunk-jpl/src/m/modeldata/interpShepherd2019.m (revision 27548) +++ /issm/trunk-jpl/src/m/modeldata/interpShepherd2019.m (revision 27548) @@ -0,0 +1,53 @@ +function dhdt=interpShepherd2019(X,Y,string,varargin) +%INTERPSHEPHERD2019 - interpolate Shepherd2019 data +% +% Available data: +% 1. dhdt_1992_1996 +% 2. dhdt_1997_2001 +% 3. dhdt_2002_2006 +% 4. dhdt_2007_2011 +% 5. dhdt_2012_2016 +% 6. dhdt_1992_2017 +% 7. uncert_1992_1996 +% 8. uncert_1997_2001 +% 9. uncert_2002_2006 +% 10. uncert_2007_2011 +% 11. uncert_2012_2016 +% 12. uncert_1992_2017 +% +% Usage: +% [dataout] = interpShepherd2019(X,Y,'dhdt_1992_2017') + +options={'dhdt_1992_1996','dhdt_1997_2001','dhdt_2002_2006','dhdt_2007_2011','dhdt_2012_2016','dhdt_1992_2017',... + 'uncert_1992_1996','uncert_1997_2001','uncert_2002_2006','uncert_2007_2011','uncert_2012_2016','uncert_1992_2017'}; +tf=strcmp(string,options); + +if ~any(tf) + disp('String not available!'); + disp('The options are:'); + disp(options); + error('String not available. See message above.'); +end + +switch oshostname(), + case {'ronne'} + nc='/home/ModelData/Antarctica/DHDTShepherd/antarctic_dhdt_5km_grid_1992_2017.nc'; + case {'recruta'} + nc='/home/santos/ModelData/CPOM_dhdt_shepherd_2019/antarctic_dhdt_5km_grid_1992_2017.nc'; + otherwise + error('machine not supported yet'); +end + +if nargin==3, + method='linear';% default +else + method=varargin{1}; +end + +xdata=double(ncread(nc,'x')); +ydata=double(ncread(nc,'y')); +data=double(ncread(nc,string))'; + +dhdt=InterpFromGrid(xdata,ydata,data,X,Y,method); + +end Index: /issm/trunk-jpl/src/m/modeldata/interpSmith2020.m =================================================================== --- /issm/trunk-jpl/src/m/modeldata/interpSmith2020.m (revision 27548) +++ /issm/trunk-jpl/src/m/modeldata/interpSmith2020.m (revision 27548) @@ -0,0 +1,51 @@ +function dhdt=interpSmith2020(X,Y,string,varargin) +%INTERPSMITH2020 - interpolate Smith2020 data +% +% Available data: +% Filtered mass-change maps, for display only, units of m(ice-equivalent)/yr: +% ais_floating_filt +% ais_grounded_filt +% gris_filt +% +% Raw mass-change maps, suitable for generation of basin-by-basin mass-change estimates, units of m(ice-equivalent)/yr: +% ais_floating +% ais_grounded +% gris +% +% Usage: +% [dataout] = interpSmith2020(X,Y,'ais_floating_filt') + +options={'ais_floating_filt','ais_grounded_filt','gris_filt','ais_floating','ais_grounded','gris'}; +tf=strcmp(string,options); + +if ~any(tf) + disp('String not available!'); + disp('The options are:'); + disp(options); + error('String not available. See message above.'); +end + +switch oshostname(), + case {'ronne'} + if strcmp(string,'gris_filt') | strcmp(string,'gris') + path='/home/ModelData/Greenland/DHDTSmith/'; + else + path='/home/ModelData/Antarctica/DHDTSmith/'; + end + case {'totten'} + if strcmp(string,'gris_filt') | strcmp(string,'gris') + path='/totten_1/ModelData/Greenland/DHDTSmith/'; + else + path='/totten_1/ModelData/Antarctica/DHDTSmith/'; + end + case {'recruta'} + path='/home/santos/ModelData/ICESat1_ICESat2_mass_change/'; + otherwise + error('machine not supported yet'); +end + +file=strcat(path,string,'.tif'); + +dhdt=interpFromGeotiff(file,X,Y); + +end