How to plot a filled contour on a map plot?

I'm working on a project where I need to plot a filled contour onto a map. Ideally, this contour plot would be partially transparent so that the map underneath can still be seen. I was previously having issues where matlab didn't like plotting the output from countourf onto a geomap, so I implemented the solution found here:

https://www.mathworks.com/matlabcentral/answers/784501-how-to-combine-geoaxes-and-contourf-plot

This, however, only gives the contour lines, and not the filled areas. I also tried using a mapplot, but I have other data that also needs to be on this figure and that didn't end up working.

Here is the function I used, getContourLineCoordinates:

https://www.mathworks.com/matlabcentral/fileexchange/74010-getcontourlinecoordinates

and here is my current code:

geoscatter(points,"y","x","filled","MarkerFaceColor","#0000ff")

for m=1:length(shapes)

geoplot(shapes(m))

end

%plotting contour map from previous calculation

contdata = contourc(xrange,yrange,Ftotnorm',100);

cTbl = getContourLineCoordinates(contdata); % from the file exchange

% Plot contour lines

hold(ax,'on')

nContours = max(cTbl.Group);

colors = autumn(nContours);

for i = 1:nContours

gidx = cTbl.Group == i;

geoplot(ax, cTbl.Y(gidx), cTbl.X(gidx),'Color', colors(i,:)); % note: x & y switched

end

% Add colorbar

colormap("autumn")

colorbar(ax)

[Dyuman Joshi]

Without the data, it is not possible to see the current output of the code nor to suggest anything.

Please attach the data you are working with, use the paperclip button to do so.

[Alan Garger]

Sorry about that, but I've attached the data. Points might have a different name than what I put in this code here.

[Star Strider]

I delete my unaccepted answers. You don’t seem to have used my solution anyway.

Answer 1

[Alan Garger]

Alright after seeing the comment from Star Strider I realized that I needed to take a different approach rather than just trying to plot the contour data straight onto the map. What I did instead was to read the points from each layer and then plot each layer as a separate geopolyshape, similar to how the solution I mentioned in the original question plotted the contour lines onto the map.

figure()
ax = geoaxes(); 
geobasemap('satellite')

%plotting contour map from calculation

%calculating contour data
contdata = contourc(xrange,yrange,Ftotnorm',100); 

Unrecognized function or variable 'xrange'.

cTbl = getContourLineCoordinates(contdata); % from the file exchange
hold(ax,'on')
nContours = max(cTbl.Group); 

%finding every unique level
ltable=cTbl.Level;
levs=unique(ltable);

gpart=[];

%for loop that splits countour data into groups and adds a NaN after every
%group, then recombines the data. NaNs are added to be able to plot shapes
%with holes
for i=1:nContours

    temp=cTbl(cTbl.Group==i,[1,3,4]);
    temp=temp{:,:};
    temp=[temp;temp(1,1),NaN,NaN];
    gpart=[gpart;temp];
    
end

%converts the data into a table for a search in the next for loop
gpart=array2table(gpart,'VariableNames',["Level","Long","Lat"]);

%creates a cell array that will be useful later
lpart=cell(length(levs),1);

%for loop that separates groups by their layer, ie, its value for the
%"Z" axis
for i=1:length(levs)
    
    temp=gpart(gpart.Level==levs(i),:);
    temp=temp{:,:};
    lpart(i)=mat2cell(temp,height(temp),3);

end

%color array that will be used to color the contour. Change "autumn" to
%change color gradient used
colors = autumn(length(lpart));

%for loop to actually plot the contour 
for i=1:length(lpart)

    %coordinates of the ith layer
    shapemat1=cell2mat(lpart(i));
    %coordinates of the next layer up, if statement for if there is no next
    %layer
    if i==length(lpart)
        shapemat2=[];
    else
        shapemat2=cell2mat(lpart(i+1));
    end
    %combines this layer with the next layer
    shapemat=[shapemat1;shapemat2];
    %plots the shape from above. By combining the ith layer with the i+1
    %layer, the ith layer get plotted and then the i+1 layer is then taken
    %out. This is why I did the NaN stuff from before, to split the data up
    %into something geopolyshape can recognize.
    cshape=geopolyshape(shapemat(:,3),shapemat(:,2));
    %actually plots the shape onto the map. alpha values are transparency
    geoplot(cshape,'FaceColor',colors(i,:),'FaceAlpha',.25,'EdgeAlpha',.5)

end

geolimits([40.4,41], [-75.9,-75]);

% Add colorbar
colormap("autumn")
colorbar(ax)

Granted, this code still depends on the use of getContourLineCoordinates which can be found here. I tried uploading cTbl so that it can output here but it didn't want to upload. Trust me though, it worked. This method is also likely very ineffiecient and unoptimized, but hey, it works. Feel free to make it better if you so please