How can I plot contour lines on the generated map?

Hi all,

Using the first section of code below, I could generate map of Alaska state. I need to draw contour lines (second section of the code) on the generated map. Can anyone help me with this regard? Please use .shp file in the zipped folder,

%%%%% generating Alaska map
S=shaperead('Location of the .shp file\tl_2018_02_anrc.shp','UseGeoCoords',true);
geoshow(S,"DisplayType","multipoint")
xlim([-179.148909,-130]);
ylim([51.214183,71.365162]);
%%%%%%%% generating contour lines based on the latitude and longitude of
%%%%%%%% weather stations within the Alaska sate
load('POFDE.mat');
data = cell2mat(POFDE);
[lat,lon] = meshgrid(unique(data(:,1)),unique(data(:,2)));
for i = 1:(size(data, 2)-2)
    figure
    I = scatteredInterpolant(data(:,[1 2]), data(:,i+2));
    contourm(lat,lon,min(1,max(0,I(lat,lon))),'ShowText','on')
    colorbar
    title(['The probability of frost depth exceedance of \Omega_{\delta} = ' num2str(i) ' feet'])
    exportgraphics(gca, sprintf('FrostPlot_%d_feet.png', i))
end

[Walter Roberson]

You do not need meshgrid() for contour() . contour() will accept vectors of coordinates.

meshgrid() or ndgrid() are useful in cases where the Z coordinate is being calculated, or the Z value is being interpolated, but these days many of the graphing functions accept marginal vectors.

Answer 1

[Voss]

In your code, you make the map then create a new figure for each contour. Those new figures aren't going to have the map; the map was in the first figure. To have a map in each figure, call geoshow after figure.

unzip('tl_2018_02_anrc.zip')

%%%%% generating Alaska map
% S=shaperead('Location of the .shp file\tl_2018_02_anrc.shp','UseGeoCoords',true);
S=shaperead('tl_2018_02_anrc.shp','UseGeoCoords',true);

%%%%%%%% generating contour lines based on the latitude and longitude of
%%%%%%%% weather stations within the Alaska sate
load('POFDE.mat');
data = cell2mat(POFDE);
[lat,lon] = meshgrid(unique(data(:,1)),unique(data(:,2)));
for i = 1:(size(data, 2)-2)
    figure
    geoshow(S,"DisplayType","multipoint")
    xlim([-179.148909,-130]);
    ylim([51.214183,71.365162]);
    I = scatteredInterpolant(data(:,[1 2]), data(:,i+2));
    contourm(lat,lon,min(1,max(0,I(lat,lon))),'ShowText','on')
    colorbar
    title(['The probability of frost depth exceedance of \Omega_{\delta} = ' num2str(i) ' feet'])
    exportgraphics(gca, sprintf('FrostPlot_%d_feet.png', i))
end

Warning: Function GEOSHOW expected DisplayType 'multipoint' to match Geometry 'polygon'. GEOSHOW is ignoring the DisplayType value.

Warning: Duplicate data points have been detected and removed - corresponding values have been averaged.

Warning: Function GEOSHOW expected DisplayType 'multipoint' to match Geometry 'polygon'. GEOSHOW is ignoring the DisplayType value.

Warning: Duplicate data points have been detected and removed - corresponding values have been averaged.

Warning: Function GEOSHOW expected DisplayType 'multipoint' to match Geometry 'polygon'. GEOSHOW is ignoring the DisplayType value.

Warning: Duplicate data points have been detected and removed - corresponding values have been averaged.

Warning: Function GEOSHOW expected DisplayType 'multipoint' to match Geometry 'polygon'. GEOSHOW is ignoring the DisplayType value.

Warning: Duplicate data points have been detected and removed - corresponding values have been averaged.

Warning: Function GEOSHOW expected DisplayType 'multipoint' to match Geometry 'polygon'. GEOSHOW is ignoring the DisplayType value.

Warning: Duplicate data points have been detected and removed - corresponding values have been averaged.

[Behrooz Daneshian]

you are right. However, scatteredInterpolant function must interpolate between the prescribed points(weather stations) within the map area. If you notice, contour lines generaly lies out of the map. I want them to be drawn within the map domain.

The attached figure is contour lines when they are not superimposed to the map.

[Walter Roberson]

The values should not be calculated outside of the shape, so the contour should rely only on data from inside the shape? Or can the values be calculated outside of the shape but you want to clip the drawing of the contour lines against the shape?

It might be a bit of a nuisance, but a shape structure can be converted into a mask, and the mask could be used to set portions of the data to NaN, so that only data from inside the shape was used to interpolate. Or go ahead and interpolate without the NaN but set NaN in the interpolated data.

At the moment I do not see a function designed to convert a geoshape to a mask, but I can see from the data structure that it would be possible with some work.

[Behrooz Daneshian]

Actually, the contours should rely on data from inside.Please look at the locations of weather stations on map on the attached figure. The probabilities(values shown by contour) is estimated at each weather stations. So contour lines must be estimatd based on the values there.

[Walter Roberson]

At approximately x = -171 y = 57 (lower left quardrant), there appears to be a single weather station. If you restrict interpolation to be from data inside the shape, then that weather station would be surrounded by a barrier of all NaN, and when you do 2D interpolation with NaN around you, the result is NaN unless you happen to query exactly at the location of the weather station. For example if you were to ask for the reading 3 cm to the east of the weather station, the interpolated result would have to be NaN because the adjacent readings would be from outside the shape, and NaN+something is NaN.