Create an area mask using a shapefile

I am trying to mask my study area using a shapefile But unable to do so. Here is the code that i was trying.. I'll be really grateful if someone can help or provide me with some alternative method to achieve it.

Thank you so much.

% Load the data
load Data.mat
shapefile = 'India_shapefile.shp';

S = shaperead(shapefile);
polygon = polyshape([S.X], [S.Y]);

% Create a logical mask
logical_mask = inpolygon(lon, lat, polygon.Vertices(:, 1), polygon.Vertices(:, 2));

% Use the logical mask to extract data from matrix
extracted_data = grid_data(logical_mask);

Answer 1

[Walter Roberson]

For polyshape objects use isinterior instead of inpolygon.

[Vidushi Koliyan]

Thanks for the reply but it didn't help me

[Walter Roberson]

You are using a shape file for India, but your lat and lon correspond mostly to the Indian Ocean up to roughly Singapore. For example Jakarta is inside your search area.

[Vidushi Koliyan]

Okay, I got it. Thank you for pointing out the mistake.

[Vidushi Koliyan]

Hello Walter, I have another question related to the previous one. While examining my netCDF files, which I used to extract my data, I noticed that my original data and the shapefile appear to be inversely mirrored rather than overlapping correctly. I can confirm that my shapefile has the correct projection. I've tried various rotation methods for my original grid, but none seem to align the two datasets properly. Can you offer any assistance or insights regarding this issue?

[Walter Roberson]

I would need something more specific -- the data file and to have you point out at least two example points for me to look at.

A common problem when examining netCDF files is that netCDF arrays are transposed compared to MATLAB conventions, and ncread does not transpose them when reading them. If the netcdf file says that the first dimension is latitude and the second dimension is longitude and you ncread() then what you get in MATLAB will have a first dimension that is longitude and a second dimension that is latitude. Another way of saying this is that arrays in netCF are indexed across rows -- netcdf_A(3,5) is third column 5th row in the file.

[Vidushi Koliyan]

Here is my dataset. It appears that in this context, latitudes are considered negative in the Northern Hemisphere. Consequently, when attempting to superimpose "borders" on the existing data, they do not align as expected. Even reversing the sign of latitudes does not resolve the issue.

[Vidushi Koliyan]

As you asked me to point out example points, if you see this data.. imagesc(lon, lat, data) hold on borders . I guess you will understand the problem. It will be really helpful if you can suggest something.

[Walter Roberson]

% Load the data
load data.mat
shapefile = fullfile('India Shapefile', 'India_shapefile.shp');

S = shaperead(shapefile);

%remove small triangles and squares
coord_lengths = arrayfun(@(s) numel(s.X), S);
not_useful_mask = coord_lengths < 10;
S(not_useful_mask) = [];

polygon = polyshape({S.X}, {S.Y});

% Create a logical mask
[Lon, Lat] = meshgrid(lon, lat);
logical_mask = reshape(isinterior(polygon, Lon(:), Lat(:)), size(Lon));

% Use the logical mask to extract data from matrix
extracted_data = test_data(logical_mask);

subplot(3,1,1);
imagesc(lon, lat, test_data); title('data'); xlabel('lon'); ylabel('lat'); axis equal; set(gca, 'YDir', 'normal');
subplot(3,1,2);
plot(polygon); title('borders'); xlabel('lon'); ylabel('lat'); axis equal;  set(gca, 'YDir', 'normal');
subplot(3,1,3);
imagesc(lon, lat, logical_mask); title('selected points'); xlabel('lon'); ylabel('lat'); axis equal; set(gca, 'YDir', 'normal');

You have to be careful with X and Y. Y is latitude, X is longitude, but you were passing in latitude (Y) into imagesc in the place it expects X coordinates. Also by default imagesc flips images, so we account for that.

Remember that extracted_data is going to be a vector, in the order that find() would report that logical_mask entries are true, which scans down columns. All of the true locations for one column will be extracted, then all of the true locations for the next column to the right, and so on. The code does not construct a "bounding box" or a convex hull. The extracted_data vector is generally not useful if you need to make any kind of spatial analysis. If you need to do something like contour the data within the bounds, then you would be better of doing something like

selected_data = test_data;

selected_data(~logical_mask) = nan;

[Vidushi Koliyan]

Thank you so much...