[27912] | 1 | function 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);
|
---|
[23873] | 7 |
|
---|
[27912] | 8 |
|
---|
[26701] | 9 | if nargin < 4
|
---|
| 10 | nanValue = 10^30;
|
---|
[27912] | 11 | fillholes = false;
|
---|
[26701] | 12 | end
|
---|
[27912] | 13 | if nargin < 5
|
---|
| 14 | fillholes = false;
|
---|
| 15 | end
|
---|
[26701] | 16 |
|
---|
[23873] | 17 | usemap = 0;
|
---|
| 18 | if license('test','map_toolbox')==0,
|
---|
| 19 | disp('WARNING: map toolbox not installed, trying house code');
|
---|
| 20 | usemap = 0;
|
---|
| 21 | elseif license('checkout','map_toolbox')==0
|
---|
| 22 | disp('WARNING: map toolbox not available (checkout failed), trying house code');
|
---|
| 23 | usemap = 0;
|
---|
| 24 | end
|
---|
| 25 |
|
---|
| 26 | if 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;
|
---|
| 33 | else
|
---|
| 34 |
|
---|
| 35 | %Get image info
|
---|
| 36 | Tinfo = imfinfo(geotiffname);
|
---|
[26701] | 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);
|
---|
[23873] | 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
|
---|
[24747] | 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
|
---|
[26701] | 71 | if nanValue > 0
|
---|
| 72 | data(find(abs(data)>=nanValue))=NaN;
|
---|
| 73 | else
|
---|
| 74 | data(find(data<=nanValue))=NaN;
|
---|
| 75 | end
|
---|
[27912] | 76 | if fillholes
|
---|
| 77 | disp('Filling holes');
|
---|
| 78 | data = inpaint_nans(data);
|
---|
| 79 | disp('done');
|
---|
| 80 | end
|
---|
[23873] | 81 | end
|
---|
| 82 |
|
---|
| 83 | dataout = InterpFromGrid(xdata,ydata,data,X,Y);
|
---|
| 84 | dataout(dataout==-9999)=NaN;
|
---|