How to subset the matrix/array with shape file?

I have matrix of 890*751 and array 890*751*5 and I want to clip/crop/subset this matrix and array with shape file. The shape file is consist on several sub-polygon inside. Each inside sub-polygon of shape file has to clip/crop/subset separately and save with unique name. E.g, we have a shape file with 4 sub-polygon (with name (a,b,c,d) and we have to get 4 clip/crop/subset(a,b,c,d) of the matrix/array with. we can read the shape file by

area = shaperead('shapefile', 'UseGeoCoords', true);
geoshow( area , 'Facecolor', 'blue')

but no idea how to crop/clip/subset the matrix/array. The sample testmatrix, testarray and shapefile are attached here for easy to understand the question. Thank you in advance for any helpe/guide

[KSSV]

YOu need to provide the lat and lons for matrices.

[Shakir Hussain]

Dear KSSV, the lat lon of matrix are attached here

[KSSV]

aaah...got it....

Answer 1

[KSSV]

Check this link for a similar question. You can follow the same code.

https://in.mathworks.com/matlabcentral/answers/403555-how-to-extract-data-through-the-shapefile-india-from-global-precipitation-data

shapefile = 'shapefile.shp' ;
S = shaperead(shapefile) ;

N = length(S) ;
for i = 1:N
    plot(S(i).X,S(i).Y)
    hold on
end
lon = load('testlon.mat') ; lon = lon.testlon ;
lat = load('testlat.mat') ; lat = lat.testlat ;

[X,Y] = meshgrid(lon,lat) ;

data = load('testarray.mat') ; data = data.testarray ;
[nx,ny,d] = size(data) ;
%%Extract data 
iwant = cell(d,N) ;
for i =1:d
    A = data(:,:,i) ;
    for j = 1:N
        idx = inpolygon(X(:),Y(:),S(i).X,S(i).Y) ;
        iwant{i,j} = A(idx) ;
    end
end

[Shakir Hussain]

Thanks your response.

shapefile = 'AoiTp.shp';
S = shaperead(shapefile);
N = length(S);
for i = 1:N
    plot(S(i).X,S(i).Y)
    hold on
end
lat = testlat;
lon = testlon;
[X,Y] = meshgrid(lon,lat);
data = testarray;
[nx,ny,d] = size(data);
%%Extract data 
iwant = cell(d,N);
for i =1:d
    A = data(:,:,i);
    for j = 1:N
        idx = inpolygon(X(:),Y(:),S(i).X,S(i).Y);
        iwant{i,j} = A(idx);
    end
end

Why I am getting the error "Index exceeds matrix dimensions". If it work, then How we can extract the data by sub-polygon of a shape file?

[Walter Roberson]

The number of entries in the shapefile, N, is not definitely the same as the size of the third dimension in testarray, especially if testarray turns out to be 4 dimensional.

[KSSV]

Which line shows error? Code worked fine with your data in my pc.

[Shakir Hussain]

I think Walter is right, the third line of code (N = length(S) is equal to 5 and iwant = cell(d,N); is extracting only 5 instead of 17.

[Shakir Hussain]

Error in ==> run at 18 idx = inpolygon(X(:),Y(:),S(i).X,S(i).Y) ;

Here is the error I am getting from the code

[KSSV]

What error?

[Shakir Hussain]

Error in ==> run at 18 idx = inpolygon(X(:),Y(:),S(i).X,S(i).Y) ; Here is the error I am getting from the code ( Index exceeds array bounds)

[KSSV]

What is i? I think i is going beyond length(S).

[Shakir Hussain]

Yes, Now it is fixed and code is working well. Last thing is how we convert the result (iwant) into 2d matrix according to the polygon extend?

[Adi Mulsandi]

Hi Shakir Hussain can you share how to convert the result (iwant) into 2d matrix according to the polygon extend?