Is there a fix for a bug in the gshhs function from the MATLAB mapping tool box? Faulty section of coast off South America.

I am creating maps of Antarctica and the surrounding Southern Ocean to show ocean chlorophyll content, so using a polar viewpoint with a stereo map projection and using the current version of MATLAB, MATLAB_R2014a. When using the gshhs function from MATLAB's mapping toolbox to create the coastlines, myself and my supervisor have found a visible fault along South America's East coast and a section of the surrounding waters. This occurs with both the course and low level resolution files, we are yet to try the higher resolutions (they obviously take longer processing times).

Has anyone else come across this before and/or does anyone have a solution to this?

An example of the maps we are producing with the specified problems is attached and the code we used to create it is given below. Any help would be much appreciated.

close all; clear all;

S = netcdf.open('S20060012006031.L3m_MO_SO_Chl_9km.Johnson_SO_Chl.nc4', 'NC_NOWRITE');

% for i = 0:10
%     netcdf.inqVar(S,i)
% end

chl = netcdf.getVar(S,0);
lat = netcdf.getVar(S,1);
lon = netcdf.getVar(S,2);
netcdf.close(S);

z = find(lon<0); lon(z) = lon(z)+360;
lon = lon';
for i = 1:720
clon(i,:) = lon;
end
for i = 1:4320
clat(:,i) = lat;
end
z = find(chl <0); chl(z) = NaN; clear z;
chl = double(chl');
% pcolor(lon, lat, log10(chl)); shading flat;

world = gshhs('gshhs_c.b');

h = figure('Color','white')
set(h,'papertype','tabloid')
set(h,'paperposition',[0 0 10 5]); %Gives exactly 2000 x 1000 pixel image
set(h,'InvertHardCopy','off'); %keep set background colors when printing to file

a = axesm('MapProjection','stereo','MapLatLimit',[-90 -30],'MapLonLimit',[0 360],'Frame','on','Grid','off','GLineStyle','-','GLineWidth',1,'MLineLocation',2,'PLineLocation',1,'MeridianLabel','on','ParallelLabel','on','PLabelLocation',[-80:10:-30],'MLabelLocation',[0:30:360],'FontSize',16,'FontWeight','Bold','LabelRotation','on')
set( gca, 'Visible', 'off' ) ;

geoshow([world.Lat], [world.Lon], 'DisplayType','polygon','FaceColor',[0.5 0.5 0.5])
hold on
pcolorm(lat, lon, log10(chl));

Thank you in advance!

Answer 1

[Chad Greene]

All of Matlab's built-in coastlines are pretty lousy, but here's a fix:

Delete world = gshhs('gshhs_c.b'); and delete the geoshow([world.Lat], [world.Lon], 'DisplayType','polygon','FaceColor',[0.5 0.5 0.5]) lines. Instead, at the very end of everything, try this:

load coast; 
patchm(lat,long,.5*[1 1 1])

I say do this at the very end because load coast loads an array named lat and we don't want to confuse it with the lat matrix you already have.

Also, two tricks that are unrelated to your question: First, I see a seam along the meridian where your chlorophyll data meets itself. Overlap that data by repeating a column of data as described here. And second, when you create your colormap, it's a nice habit to use jet(256) instead of the default jet, which only creates 64 colors. The difference is an aesthetic preference of mine and is not integral to the data display, but the smooth color gradient helps bring Matlab graphics into the 21st century.

[Chad Greene]

Also, setting up those polar stereographic axes can be a pain. For a simpler solution, I recommend ncpolarm.

[Amy]

Thank you so much Chad!

My revised code is pasted below, with the figure it produces attached. My next issue that I'm hoping you may be able to help with is how to plot the East and West coordinate labels so that they are just outside of edge of the map, rather than cutting through the line of latitude?

    close all; clear all;

S = netcdf.open('S20060012006031.L3m_MO_SO_Chl_9km.Johnson_SO_Chl.nc4', 'NC_NOWRITE');

% for i = 0:10
%     netcdf.inqVar(S,i)
% end

chl = netcdf.getVar(S,0);
chllat = netcdf.getVar(S,1);
chllon = netcdf.getVar(S,2);
netcdf.close(S);

z = find(chllon<0); chllon(z) = chllon(z)+360;
chllon = chllon';
% for i = 1:720
% clon(i,:) = lon;
% end
% for i = 1:4320
% clat(:,i) = lat;
% end
z = find(chl <0); chl(z) = NaN; clear z;
chl = double(chl');
% pcolor(lon, lat, log10(chl)); shading flat;

h = figure('Color','white')
set(h,'papertype','tabloid')
set(h,'paperposition',[0 0 10 5]); %Gives exactly 2000 x 1000 pixel image
set(h,'InvertHardCopy','off'); %keep set background colors when printing to file

a = axesm('MapProjection','stereo','MapLatLimit',[-90 -30],'MapLonLimit',[0 360],'Frame','on','Grid','on','GLineStyle','-','GLineWidth',1,'MLineLocation',30,'PLineLocation',20,'MeridianLabel','on','ParallelLabel','on','PLabelLocation',[-80:20:-40],'MLabelLocation',[-180:30:150],'FontSize',14,'FontWeight','normal','LabelRotation','on')
set( gca, 'Visible', 'off' ) ;

hold on
pcolorm(chllat, chllon, log10(chl));
hold on 
load coast; 
coastlat = lat; clear lat; coastlon = long; clear long;
patchm(coastlat,coastlon,0.5*[1 1 1])

% shading interp

set(gca, 'clim', [log10(0.1) log10(10)]); % sets axes limits for chl

cb = colorbar('horiz');

set(cb, 'xtick',[log10(0.1) log10(0.3) log10(0.5) log10(1) log10(5) log10(10)],'xticklabel',[0.1 0.3 0.5 1 5 10], 'fontsize', 16,'fontweight', 'bold');%for chl
set(h, 'PaperPositionMode', 'auto');
print -djpeg -r600 test.jpgif

Any help is much appreciated! Thanks, Amy

[Chad Greene]

Follow it all with:

 setm(a,'mlabelparallel',-35)

to put the meridian labels at 35 degrees south. Unfortunately if you try to put them farther north, they'll all scrunch up at the equator. In fact, the reason the outside edge of the map looks like a bold line is because all the parallels are stacked right there.

[Stephen23]

"I say do this at the very end because load coast loads an array named lat and we don't want to confuse it with the lat matrix you already have"

A much simpler and more reliable solution is to always load into a structure (which is always recommended anyway), then you will never have this problem at all:

S = load('coast'); 
patchm(S.lat,S.long,.5*[1 1 1])

Answer 2

[Kelly Kearney]

Just to clarify, what you're seeing isn't a bug in gshhs so much as a bug in geoshow, patchm, patchesm, etc.; all of Matlab's functions designed to plot polygons on maps seem to fall apart for all but the simplest of polygons. The coast.mat polygons generally pass the test, but no other coastline dataset does so reliably, in my experience.

My usual workaround is to triangulate any patches that I'm planning to plot to a map axis, as I explain in this thread. I've submitted bug reports and enhancement requests related to this over the past 8 or so years, but plotting maps in Matlab still seems to be a bit more of an art than a science.

-Kelly