source: issm/trunk-jpl/src/m/contrib/morlighem/modeldata/interpFromGeotiff.m@ 26701

Last change on this file since 26701 was 26701, checked in by Cheng Gong, 3 years ago

ADD:interpret MEaSUREs Greenland Ice velocity to mesh

File size: 1.8 KB
Line 
1function dataout = interpFromGeotiff(geotiffname,X,Y,nanValue),
2
3if nargin < 4
4 nanValue = 10^30;
5end
6
7usemap = 0;
8if license('test','map_toolbox')==0,
9 disp('WARNING: map toolbox not installed, trying house code');
10 usemap = 0;
11elseif license('checkout','map_toolbox')==0
12 disp('WARNING: map toolbox not available (checkout failed), trying house code');
13 usemap = 0;
14end
15
16if usemap,
17 [data,R] = geotiffread(geotiffname);
18 data=double(flipud(data));
19 xdata=R.XLimWorld(1):R.DeltaX:R.XLimWorld(2); xdata=xdata(:);
20 xdata =(xdata(1:end-1)+xdata(2:end))/2;
21 ydata=R.YLimWorld(2):R.DeltaY:R.YLimWorld(1); ydata=flipud(ydata(:));
22 ydata =(ydata(1:end-1)+ydata(2:end))/2;
23else
24
25 %Get image info
26 Tinfo = imfinfo(geotiffname);
27 N = Tinfo(1).Width;
28 M = Tinfo(1).Height;
29 dx = Tinfo(1).ModelPixelScaleTag(1);
30 dy = Tinfo(1).ModelPixelScaleTag(2);
31 minx = Tinfo(1).ModelTiepointTag(4);
32 maxy = Tinfo(1).ModelTiepointTag(5);
33
34 %Generate vectors
35 xdata = minx + dx/2 + ((0:N-1).*dx);
36 ydata = maxy - dy/2 - ((M -1:-1:0).*dy);
37
38 %Read image
39 if 1
40 assert(dx>0); assert(dy>0);
41 ydata = fliplr(ydata);
42
43 %Get pixels we are interested in
44 offset=2;
45 xmin=min(X(:)); xmax=max(X(:));
46 posx=find(xdata<=xmax);
47 id1x=max(1,find(xdata>=xmin,1)-offset);
48 id2x=min(numel(xdata),posx(end)+offset);
49
50 ymin=min(Y(:)); ymax=max(Y(:));
51 posy=find(ydata>=ymin);
52 id1y=max(1,find(ydata<=ymax,1)-offset);
53 id2y=min(numel(ydata),posy(end)+offset);
54
55 data = double(imread(geotiffname,'PixelRegion',{[id1y,id2y],[id1x,id2x]}));
56 xdata=xdata(id1x:id2x);
57 ydata=ydata(id1y:id2y);
58 else
59 data=double(flipud(imread(geotiffname)));
60 end
61 if nanValue > 0
62 data(find(abs(data)>=nanValue))=NaN;
63 else
64 data(find(data<=nanValue))=NaN;
65 end
66
67end
68
69dataout = InterpFromGrid(xdata,ydata,data,X,Y);
70dataout(dataout==-9999)=NaN;
Note: See TracBrowser for help on using the repository browser.