Taking subset of geotiff values based on longitude and latitude.

Hi all

I have been trying to pull a subset of values from a series of .tif rasters based on longitude and latitude bounds. The rasters in questions are at a global scale where the data of interest lies between the degrees of 49 - 72 long , 31 - 51 lat. (Which are stored in variables X & Y in the following code)

The rasters in question have a WGS 84 datum and a pixel scale of 2.5' x 2.5'.

I have found a solution which is inelegant so was hoping that someone would be able to be of assistance.

   %%Open and Index files containing Raster images
  tifs=dir('C:\Users\c3083\OneDrive - The University Of Newcastle\FYP\Nathan DICK\MATLAB Model\RainfallM_2.5');
  Srch='RainRast';
  Rnm=extractfield(tifs,'name');
  tmp=find(contains(Rnm,Srch));
  Rnm=Rnm(tmp);
  Rnm=natsortfiles(Rnm); % function dl from web to sort name strings
  stps=numel(Rnm); 

%%lat long as matrix references
LatR=([min(Y) max(Y)])/r.CellExtentInLatitude;
LongR=([min(X) max(X)]+180)/r.CellExtentInLongitude;

for i = 1:stps
    [Pt, r]=geotiffread(char(Rnm(i)));
    P(:,:,i)=cat(3,Pt(LatR(1):LatR(2),LongR(1):LongR(2)));
end

view=P(:,:,1);

clear Pt tifs Rnm Srch 

%%Junk code which hasn't worked but shows posting wasn't a first resort .
% rr=geotiffinfo('RainRast_01.tif');
% [Pt, r]=geotiffread(char(Rnm(1)));
% [xx, yy]=pixcenters(rr);
% [mshX, mshY]=meshgrid(xx,yy);
% mask=inpolygon(mshX, mshY, X, Y);
%  P=zeros(LatR(2)-LatR(1)+1,LongR(2)-LongR(1)+1,stps);
% % [Pt, r]=geotiffread(char(Rnm(i)));
% % [RasX, RasY]=pix2map(info.r, rows
%vect=combvec(xx,yy);
%[mshX, mshY]=meshgrid(xx,yy);
% Pt=union(Pt(mshX>=min(X)& mshY<=max(X),mshY>=min(Y)& mshX<=max(Y)));
% [mshX,~, destcol] = find(xx>=min(X)& xx<=max(X)); %Set col as latitude values
% [mshY, ~, destrow] = unique(Loc(:, 1)); %set row as Longitudinal Values

% [C,xx,yy]=union(Pt);
% 
% surf(mshX, mshY, Pt)
% axis([-inf inf -inf inf 0 inf])
%P(:,:,1)=cat(3,find(Pt(mshY>=min(X)& mshY<=max(X),mshX>=min(Y)& mshX<=max(Y))));

% [mshX,~, destcol]=Pt(xx>=min(X)& xx<=max(X));
% [mshY, ~, destrow] = y1(:); %set row as Longitudinal Values

Answer 1

[KSSV]

Create your area of interest/ subset grid using meshgrid fom the bounds of lat and lon.. then use interp2 to get the respective Z data.

[Emily T. Griffiths]

Could you please elaborate on this?

%Subsetting grid
[Xq,Yq]=meshgrid(-5.2:0.05:13.2,50.9:0.025:62.1); 
%Input larger data.
[data, R] = geotiffread(geoTIFF_filename);
info=geotiffinfo(geoTIFF_filename);
Vq = interp2(?,?,data,Xq,Yq)  %<-- What goes here?

Is it R.LongitudeLimits, etc? I tried that and it doesn't work:

Error using  .' 
Transpose on ND array is not defined. Use PERMUTE instead.

Error in interp2 (line 132)
        V = V.';

Permute is not what I want to do. I feel like this is a very simple thing that I am having a very hard time figuring out. I have Matlab 2017b, so I don't have most of the newly available functions that appear to do this seemlessly.

[KSSV]

R has a range of lat and lon i.e the bounding box. Create your main grid with that extents and the size of data.