Contour plot around edge on map

4 views (last 30 days)
Hi
I have some data I have to plot on a world map. It has to be a contour.
This is what i have tried:
lat = data(:, 1);
lon = data(:, 2);
z = data(:, 3);
F = scatteredInterpolant(lat, lon, z);
[lat, lon] = meshgrid(linspace(-90, 90, 200), linspace(-180, 180, 200));
z = F(lat, lon);
figure;
worldmap world;
load coastlines;
plotm(coastlat, coastlon);
hold on;
contourm(lat, lon, z);
colorbar;
This do give a world map with a contour on it. But the edges are a problematic area. It seems to me that the intpolation cannot "go around corners", which makes sense. I am not sure if the contourm function can go around the edge if the correct interpolation is made.
Do anyone know how to do a better interpolation on a world map?
  3 Comments
Michael Madelaire
Michael Madelaire on 4 Oct 2017
Thanks for the hint.
I have tried a couple of different things now: Changing the data from lat/lon to a regular data grid with geo Raster. And defining worldmap with my data.
[Z, refvec] = geoloc2grid(lat, lon, z, 0.5);
figure;
worldmap(Z, refvec);
plotm(coastlat, coastlon, 'k', 'LineWidth', 1.2);
hold on;
contourm(Z, refvec);
This is not working and I am out of ideas. Got any other ideas?
Michael Madelaire
Michael Madelaire on 4 Oct 2017
I found an "okay" solution. As I suspected, in the main question, the interpolation was to blame. If the data goes from [-180, 180] or [0, 360] the interpolation will not create a discontinuity on the map between (-180, 180) or (0, 360). In order to force the interpolation to take the "other side" into account I extended the data:
lat = [lat; lat; lat];
lon = [lon-360; lon, lon+360] % Easiest to do if your longitude is [0,360].
value = [value; value; value];
% The Interpolation will now take the other side into account.
% This is not a pretty solution and will create noise in the data.
% But one can argue that the repeated datapoints are so fare away
% that their effect is minimal to none.
F = scatteredInterpolant(lat, lon, value, 'natural');
% The meshgrid function will then "cut out" original data.
[lat, lon] = meshgrid(linspace(-90, 90, 500), linspace(0, 360, 500));
Z = F(lat, lon);
figure;
worldmap world;
plotm(coastlat, coastlon, 'k', 'LineWidth', 1.2);
hold on;
contourm(lat, lon, Z);

Sign in to comment.

Accepted Answer

Michael Madelaire
Michael Madelaire on 24 Oct 2017
I found an "okay" solution. As I suspected, in the main question, the interpolation was to blame. If the data goes from [-180, 180] or [0, 360] the interpolation will not create a discontinuity on the map between (-180, 180) or (0, 360). In order to force the interpolation to take the "other side" into account I extended the data:
lat = [lat; lat; lat];
lon = [lon-360; lon, lon+360] % Easiest to do if your longitude is [0,360].
value = [value; value; value];
% The Interpolation will now take the other side into account.
% This is not a pretty solution and will create noise in the data.
% But one can argue that the repeated datapoints are so fare away
% that their effect is minimal to none.
F = scatteredInterpolant(lat, lon, value, 'natural');
% The meshgrid function will then "cut out" original data.
[lat, lon] = meshgrid(linspace(-90, 90, 500), linspace(0, 360, 500));
Z = F(lat, lon);
figure;
worldmap world;
plotm(coastlat, coastlon, 'k', 'LineWidth', 1.2);
hold on;
contourm(lat, lon, Z);

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!