It seems like you are asking for a way to identify and retrieve every instance where the temperature value is exactly 25 within the 4D array and to return the lon, lat, depth, and years values that correspond to those locations within the array. If my interpretation of the question is correct, you're not looking to reshape or reduce the dimensions of the array (as squeeze might imply).
0. Create 4D demo array that includes values of 25 and define the 4 dimensions
rng default
A = randi(25,5,4,6,3);
lat = [0 45 90 -45 -90]; % dim 1
lon = [50 100 150 180]; % dim 2
depth = [10 20 30 40 50 60]; % dim 3
year = [1900 1950 2000]; % dim 4
criticalValue = 25;
linIdx = find(A==criticalValue);
2. Convert the linear index to subscript indices (row, column, 3rd-dim and 4th-dim indices)_
[row,col,dim3Idx,dim4idx] = ind2sub(size(A),linIdx);
3. Return the lat, lon, depth, and year for each value of 25 in the array. This assumes the 1st dimension is latitude, 2nd is longitude, 3rd is depth, and 4th is year.
lats = lat(row)
lons = lon(col)
depths = depth(dim3Idx)
years = year(dim4idx)
If the goal is to remove values of 25 from the 4D array, there are a couple of options. Assuming the values are spead throughout the array, you can't remove the data without changing the shape and size of the array and that will make it difficult to understand the lat, lon, depth, and year definitions for each dimension. Instead, you could either replace those values with missing a value indicator (NaN) or you can use an interpolation techique to replace those value. To replace with NaNs, do step 1 above to compute linIdx and then A(linIdx)=NaN;.
