How can I remove contours outside the US border?

Hello all,

Using the code below, I am going to draw contour lines showing the probability that frost depth exceeds 1 foot accros the US. Using your guidance, I used masking method in order to remove contour lines outside the US border. Althogh I could improve my contour map to some extent, there are still some contour lines lying outside the US map. Can anyone help me with this regard?

Compelelemtary explainations can be seen on comments in the code.

clear
close all
clc
%%%% Loading longitude(X) and latitude(Y) vectors of US border points which
%%%% are captured by ArcGIS
load("US_boundary_X.mat");
load("US_boundary_Y.mat");
%%%Prob1 is a cell array that its first column, second column, and third
%%%column are latitude, longitude of weather stations, and probability that
%%%frost depth exceeds 1 foot, respectively
load("Prob1.mat");
%%%%%determining bounding box for the US 
State_name=shaperead('usastatehi', 'UseGeoCoords',true);
Bounding_box={State_name.BoundingBox};
Coordinate_table=zeros(size(Bounding_box,2),4);
for i=1:size(Bounding_box,2)
    Coordinate_table(i,:)=[Bounding_box{i}(:,1)',Bounding_box{i}(:,2)'];
end
min_lon=min(Coordinate_table(:,1));
max_lon=max(Coordinate_table(:,2));
min_lat=min(Coordinate_table(:,3));
max_lat=max(Coordinate_table(:,4));
%%%%%%%%creating mesh based on the bounding box
[lon,lat]=meshgrid(linspace(min_lon,max_lon,50),linspace(min_lat,max_lat,50));
%%%%creating a polygon based on the US border coordinates(X,Y) and determining
%%%%which nodes of the mesh will be located inside the US polygon 
pgon = polyshape(X, Y);
idx = isinterior(pgon,lon(:),lat(:));
%%%%drawing a raw map of the US so that we can plot contours on that
data = cell2mat(Prob1);
figure
usamap("conus");
states = readgeotable("usastatelo.shp");
rows = states.Name ~= "Alaska" & states.Name ~= "Hawaii";
states = states(rows,:);
h = height(states);
faceColors = makesymbolspec("Polygon",...
    {'INDEX',[1 h],'Facecolor','#FFFFFF'}); 
geoshow(states,"DisplayType","polygon","SymbolSpec",faceColors)
%%%% Interpolaing using scatteredInterpolant function. We know the
%%%% probability that frost depth exceeds 1 foot in the locations of weather
%%%% stations. Now, to draw contour lines, we need to use interpolation to
%%%% estimate the probability of frost depth exceedance of 1 foot in all
%%%% nodes of the generated mesh
I = scatteredInterpolant(data(:,[2 1]), data(:,3));
z1=I(lon(:),lat(:));       
z1(~idx) = nan;
z1=reshape(z1,size(lon));
%%%% draw contours
contourm(lat,lon,min(1,max(0,z1)));
caxis([0, 1]) 
colorbar

[Walter Roberson]

See also earlier question from the same poster in which the target was specifically the dry land borders

https://www.mathworks.com/matlabcentral/answers/1950638-how-can-i-get-us-dry-land-border-coordinates-lat-and-lon?s_tid=srchtitle