Selected Land class data plot

I wanted to make the land class plot

for four selected datasets.

load('Land_Class_by_Percentage_3646N_8269W.mat')
[lonmat,latmat]=meshgrid(lonv_landclass,latv_landclass);
iwant = find(lonmat <-82.05 | lonmat >-69 | latmat >46.05 | latmat <36);
iLand_Class(iwant) = [];

iA = reshape(iLand_Class, 201,261);
lonv_landclass(lonv_landclass<-82.05 | lonv_landclass>-69) = [];
latv_landclass(latv_landclass>46.05 | latv_landclass<36.03) = [];
[lonmat,latmat]=meshgrid(lonv_landclass,latv_landclass);

%%
latlim=[36 46];
lonlim=[-82 -69];
% figure to see
figure
ax = worldmap(latlim, lonlim);
states = shaperead('usastatehi','UseGeoCoords', true, 'BoundingBox', [lonlim', latlim']);
geoshow(ax, states, 'FaceColor', [1 1 1]);
surfacem(latmat,lonmat,iA);
title('Land Classification')
colormap;
map2 = colormap; map2( 1, : ) = 1; colormap(map2);
colorbar

iwant2 = find(iA ~=5 & iA ~=13 & iA ~= 14& iA ~= 15);
iA2 =  (iA);
iA2(iwant2) = NaN;
iA2 =(reshape(iA2, size(iA)));

iA2(iA2==5) = 1;
iA2(iA2==13) = 2;
iA2(iA2==14) = 3;
iA2(iA2==15) = 4;
% figure;
% imagesc(iA2)
% get current colormap
map = colormap;
%  adjust for number of colors you want
rows = uint16(linspace(1, size(map,1), 4)) ;
map = map(rows, :);

%  apply   colormap
colormap(map);

labels = {'broadleaf forest','croplands','urban', 'vegetation'};
lcolorbar(labels,'fontweight','bold');

%%
latlim=[36 46];
lonlim=[-82 -69];
% figure to see
figure
ax = worldmap(latlim, lonlim);
states = shaperead('usastatehi','UseGeoCoords', true, 'BoundingBox', [lonlim', latlim']);
geoshow(ax, states, 'FaceColor', [1 1 1]);
surfacem(latmat,lonmat,iA2);
title('Land Classification')
%# get current colormap
map = colormap;

%# adjust for number of colors you want
rows = uint16(linspace(1, size(map,1), 4)) ;
map = map(rows, :);
%# and apply the new colormap
colormap(map);
% colorbar
% caxis([1 4])
labels = {'broadleaf forest','croplands','urban', 'vegetation'};
lcolorbar(labels,'fontweight','bold');

%%
% test to display the land classes
for jj = 1:length(classification)
   
    disp([num2str(jj), ':', classification{jj}])
end