Index: /issm/trunk-jpl/src/m/contrib/chenggong/interpFromMEaSUREsGeotiff.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/chenggong/interpFromMEaSUREsGeotiff.m	(revision 26701)
+++ /issm/trunk-jpl/src/m/contrib/chenggong/interpFromMEaSUREsGeotiff.m	(revision 26701)
@@ -0,0 +1,40 @@
+% This function calls src/m/contrib/morlighem/modeldata/interpFromGeotiff.m for multiple times to load all avaliable 
+% tif data in  /totten_1/ModelData/Greenland/VelMEaSUREs/Jakobshavn_2008_2021/ within the given time period (in decimal years)
+% For some reason, each .tif file in this folder contains two sets of data, only the first dataset is useful
+
+function dataout = interpFromMEaSUREsGeotiff(X,Y,Tstart,Tend)
+
+foldername = '/totten_1/ModelData/Greenland/VelMEaSUREs/Jakobshavn_2008_2021/';
+
+% get the time info from file names
+templist = dir([foldername,'*.meta']);
+Ndata = length(templist);
+dataTstart = zeros(Ndata,1);
+dataTend = zeros(Ndata,1);
+
+for i = 1:Ndata
+	tempConv = split(templist(i).name, '_');
+	% follow the naming convention
+	dataPrefix(i) = join(tempConv(1:5), '_');
+	dataTstart(i) = date2decyear(datenum(tempConv{3}));
+	dataTend(i) = date2decyear(datenum(tempConv{4}));
+end
+disp(['  Found ', num2str(Ndata), ' records in ', foldername]);
+disp(['    from ', datestr(decyear2date(min(dataTstart)),'yyyy-mm-dd'), ' to ', datestr(decyear2date(max(dataTend)),'yyyy-mm-dd') ]);
+
+
+% find all the data files with Tstart<=t<=Tend
+dataInd = (dataTend>=Tstart) & (dataTstart<=Tend);
+disp([' For the selected period: ', datestr(decyear2date((Tstart)),'yyyy-mm-dd'), ' to ', datestr(decyear2date((Tend)),'yyyy-mm-dd'), ', there are ', num2str(sum(dataInd)), ' records' ]);
+
+dataToLoad = dataPrefix(dataInd);
+TstartToload = dataTstart(dataInd);
+TendToload = dataTend(dataInd);
+
+for i = 1:length(dataToLoad)
+	dataout(i).vx = interpFromGeotiff([foldername, dataToLoad{i}, '_vx_v04.0.tif'], X, Y, 2e9);
+	dataout(i).vy = interpFromGeotiff([foldername, dataToLoad{i}, '_vy_v04.0.tif'], X, Y, 2e9);
+	dataout(i).vel = interpFromGeotiff([foldername, dataToLoad{i}, '_vv_v04.0.tif'], X, Y, -1);
+	dataout(i).Tstart = TstartToload(i);
+	dataout(i).Tend = TendToload(i);
+end
Index: /issm/trunk-jpl/src/m/contrib/morlighem/modeldata/interpFromGeotiff.m
===================================================================
--- /issm/trunk-jpl/src/m/contrib/morlighem/modeldata/interpFromGeotiff.m	(revision 26700)
+++ /issm/trunk-jpl/src/m/contrib/morlighem/modeldata/interpFromGeotiff.m	(revision 26701)
@@ -1,3 +1,7 @@
-function dataout = interpFromGeotiff(geotiffname,X,Y),
+function dataout = interpFromGeotiff(geotiffname,X,Y,nanValue),
+
+if nargin < 4
+	nanValue = 10^30;
+end
 
 usemap = 0;
@@ -21,10 +25,10 @@
 	%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);
+	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
@@ -55,5 +59,9 @@
 		data=double(flipud(imread(geotiffname)));
 	end
-	data(find(abs(data)>10^30))=NaN;
+	if nanValue > 0
+		data(find(abs(data)>=nanValue))=NaN;
+	else 
+		data(find(data<=nanValue))=NaN;
+	end
 
 end