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

5 views (last 30 days)
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!

Accepted Answer

Chad Greene
Chad Greene on 22 Jul 2014
Edited: Chad Greene on 22 Jul 2014
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.
  4 Comments
Chad Greene
Chad Greene on 28 Jul 2014
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
Stephen23 on 18 May 2017
Edited: Stephen23 on 18 May 2017
"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])

Sign in to comment.

More Answers (1)

Kelly Kearney
Kelly Kearney on 28 Jul 2014
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

Products

Community Treasure Hunt

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

Start Hunting!