1 | function dataout = interpFromGeotiff(geotiffname,X,Y,nanValue),
|
---|
2 |
|
---|
3 | if nargin < 4
|
---|
4 | nanValue = 10^30;
|
---|
5 | end
|
---|
6 |
|
---|
7 | usemap = 0;
|
---|
8 | if license('test','map_toolbox')==0,
|
---|
9 | disp('WARNING: map toolbox not installed, trying house code');
|
---|
10 | usemap = 0;
|
---|
11 | elseif license('checkout','map_toolbox')==0
|
---|
12 | disp('WARNING: map toolbox not available (checkout failed), trying house code');
|
---|
13 | usemap = 0;
|
---|
14 | end
|
---|
15 |
|
---|
16 | if 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;
|
---|
23 | else
|
---|
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 |
|
---|
67 | end
|
---|
68 |
|
---|
69 | dataout = InterpFromGrid(xdata,ydata,data,X,Y);
|
---|
70 | dataout(dataout==-9999)=NaN;
|
---|