source: issm/trunk-jpl/src/m/modeldata/interpFromGeotiff.m@ 27912

Last change on this file since 27912 was 27912, checked in by Mathieu Morlighem, 18 months ago

CHG: added fill holes option

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