You are trying to perform what is called a Zonal Statistics in the "GIS world". Do you really need to do it with MATLAB, or could you perform it as a one time operation in e.g. ArcGIS?
If you need to do it in MATLAB, one way to go is to define a set of compartment IDs for your data points, and to use ACCUMARRAY to get whatever statistics you need by compartment. This is comparable to what is explained here. Compartment IDs can be computed based on coordinates, e.g
compartmentID = colID + (rowID-1) * nRows ;
where compartmentID, rowID, and colID have the same size as your data set (for one year), rowID is computed based on latitudes, and colID is computed based on longitudes.
Once you have compartment IDs (1,2,..), you can get a statistics as follows
data_yr = data(:,:,1) ; % Example for year 1.
compMean = accumarray( compartmentID(:), data_yr(:), [], @mean ) ;
Note that it is sometimes more efficient to compute a mean by dividing sums..
compSum = accumarray( compartmentID(:), data_yr(:) ) ;
compCount = accumarray( compartmentID(:), ones( numel( compartmentID ), 1 )) ;
compMean = compSum ./ compCount ;
=== Illustration =====================
We use a 4 by 6 data set to experiment with:
>> data_yr = 20 + randi( 10, 4, 6 )
data_yr =
27 27 23 27 25 22
28 22 21 24 24 25
28 28 21 30 28 25
24 21 29 21 28 27
Assume that the first dimension corresponds to latitudes (1 to 4) and that the second dimension corresponds to longitudes (1 to 6). We want to get a statistics on 2 by 2 blocks/compartments, so the upper left
27 27
28 22
are in compartment 1 and so on. If we had the following array (which defines the compartment ID for each data point) available
>> compartmentID
compartmentID =
1 1 2 2 3 3
1 1 2 2 3 3
4 4 5 5 6 6
4 4 5 5 6 6
we could easily get a zonal statistics using ACCUMARRAY (I let you read the doc and experiment a bit with this function):
>> compMean = accumarray( compartmentID(:), data_yr(:), [], @mean )
compMean =
26.0000
23.7500
24.0000
25.2500
25.2500
27.0000
and you can check e.g. that the first element is the mean for compartment 1. Now the only thing left is to find a way to build compartmentID. One way to achieve this is as follows (and again, I let you experiment with it to understand):
>> rowId = ceil( (1 : size(data_yr, 1)) / 2 )
rowId =
1 1 2 2
where denom. 2 is the block size along dim 1.
>> colId = ceil( (1 : size(data_yr, 2)) / 2 )
colId =
1 1 2 2 3 3
where denom. 2 is the block size along dim 2.
>> [colID, rowID] = meshgrid( colId, rowId )
colID =
1 1 2 2 3 3
1 1 2 2 3 3
1 1 2 2 3 3
1 1 2 2 3 3
rowID =
1 1 1 1 1 1
1 1 1 1 1 1
2 2 2 2 2 2
2 2 2 2 2 2
And finally
>> compartmentID = colID + (rowID - 1) * max(colId)
compartmentID =
1 1 2 2 3 3
1 1 2 2 3 3
4 4 5 5 6 6
4 4 5 5 6 6
