% Longitude limits of plot
% Can set lower bound to negative value for a plot that
% spans the prime meridian
lons = [120, 180];

% Latitude limits of plot
lats = [20, 50];

% Density of wind vectors
% must be a positive integer
% smaller = more dense
skip = 1;

% Scale of wind vectors
% larger = longer
scale = 1.5;

% Load near-surface currents
lat = double(ncread('world_oscar_vel_clim_1993_2012.nc', 'latitude'));
lon = double(ncread('world_oscar_vel_clim_1993_2012.nc', 'longitude'));
u = double(ncread('world_oscar_vel_clim_1993_2012.nc', 'u'));
v = double(ncread('world_oscar_vel_clim_1993_2012.nc', 'v'));
umag = sqrt(u.^2 + v.^2);

% Set up map
figure('Position', [0, 0, 800, 400])
h = worldmap(lats, lons);
axis off;
setm(h, 'MapProjection','pcarree', 'MLabelParallel', 'south', ...
    'MLineLocation', 2, 'PLineLocation', 2, 'GLineStyle', '-');

% Plot sea surface height
contourfm(lat, lon, 1e2*umag', 0:10:100);
colormap(jet);
caxis([0, 100])
h = colorbar('northoutside');
ylabel(h, 'Surface current speed (cm/s)')

% Plot coastlines
load coastlines;
plotm(coastlat, coastlon, 'k');
tightmap;

% Save figure
exportgraphics(gca,'surface_current_speed.png')