How to add coastlines to plots

Hi,

Basically i am plotting Ekman velocities in the Arctic region from wind speed data stored in a netCDF file with variables lon, lat, and the wind speed data. First i plot the global wind speed, then calculate the Ekman velocity and plot global Ekman velocity but limiting the y axis (the latitude) so the resulting pcolor plot only focuses on the Arctic region. I am now trying to add a coastline to this plot of Ekman velocity in the Arctic region, i've tried many different toolboxes now like the m_map toolbox and the Arctic coastlines toolbox but none of it is working, mainly due to the fact that i don't really know how to use these toolboxes properly and i can't find much help anywhere online.

If anyone knows any simple solutions to adding a coastline to a plot, or can show me how to use these toolboxes it would be greatly appreciated. My Matlab code is as follows:

if true
clear all;
close all;
%%constants
P=1029; % sea water density (kgm^-3)
Pa=1.22; % air density (kgm^-3)
Cd=0.0016; % drag coefficient taken from Petty et al., 2017
%%U wind component
Uncid = netcdf.open('U_surface_wind_global.nc'); %Monthly means surface U-wind entire Ocean NCEP/NCAR reanalysis
Ulat=netcdf.getVar(Uncid,0);
whos Ulat;
Ulon = netcdf.getVar(Uncid,1);
whos Ulon;
Utime = netcdf.getVar(Uncid,2);
whos Utime;
Uwind = netcdf.getVar(Uncid,3);
whos Uwind;
Uwind = permute(Uwind, [2 1 3]);
figure
pcolor(Ulon', Ulat, double(Uwind(:,:,1))); shading interp; %plotting global U-wind speed
c=colorbar;
set(get(c,'title'),'string','m/s','fontsize',20);
%title('Global surface U-wind','fontsize',12);
xlabel('Longitude','fontsize',20);
ylabel('Latitude','fontsize',20);
%ylim ([60 85]);
%%V wind component
Vncid = netcdf.open('V_surface_wind_global.nc'); %Monthly means surface V-wind entire Ocean NCEP/NCAR reanalysis
Vlat=netcdf.getVar(Vncid,0);
whos Vlat;
Vlon = netcdf.getVar(Vncid,1);
whos Vlon;
Vtime = netcdf.getVar(Vncid,2);
whos Vtime;
Vwind = netcdf.getVar(Vncid,3);
whos Vwind;
Vwind = permute(Vwind, [2 1 3]);
figure
pcolor(Vlon', Vlat, double(Vwind(:,:,1))); shading interp; %plotting global V-wind speed
c=colorbar;
set(get(c,'title'),'string','m/s','fontsize',20);
%title('Global surface V-wind','fontsize',12);
xlabel('Longitude','fontsize',20);
ylabel('Latitude','fontsize',20);
%ylim ([60 85]);
%%Calculating coriolis parameter f 
f = sw_f(Ulat);  
f= repmat(f, [1, 144, 841]);
%%calcualting distances for dy and dx 
% using sw_dist package which calculates distance between two latitude and longitude coordinates
% allowing gradients to be calculated
dy =  sw_dist([0 2.5], [0 0],'km')*1000;
doublelat = repmat(Ulat, 1, 2)'; 
dx =  sw_dist(doublelat(:),repmat([0; 2.5], 73,1),'km')*1000;
dx = interp1(1:72, dx(2:2:end), 0.5:72.5)';
%%calculating gradient of Uwind
[Ux,Uy]=gradient(Uwind);
Ux = Ux./repmat(dx, [1, 144, 841]);
Uy = Uy./dy;
figure
pcolor(Ulon', Ulat, double(Ux(:, :, 1))); shading interp; %plotting gradients of global Uwind
c=colorbar;
set(get(c,'title'),'string','m/s','fontsize',20);
%title('Global surface U-wind gradient','fontsize',12);
xlabel('Longitude','fontsize',20);
ylabel('Latitude','fontsize',20);
%ylim ([70 90]); %focusing (zooming in) on Arctic region ie. latitude above
%70deg
%%Calculating gradient of Vwind
[Vx,Vy]=gradient(Vwind);
Vx = Vx./repmat(dx, [1, 144, 841]);
Vy = Vy./dy;
figure
pcolor(Vlon', Vlat, double(Vy(:, :, 1))); shading interp; %plotting gradients of global Vwind
c=colorbar;
set(get(c,'title'),'string','m/s','fontsize',20);
%title('Global surface V-wind gradient','fontsize',12);
xlabel('Longitude','fontsize',20);
ylabel('Latitude','fontsize',20);
%ylim ([70 90]); %focusing (zooming in) on Arctic region ie. latitude above
%70deg
%%Calculating global Uwind stress
Usq = Uwind.*abs(Uwind); % Uwind speed squared
Tu=Pa*Cd*Usq; % Tu=global Uwind stress
Tu=permute(Tu,[1 2 3]);
figure
pcolor(Ulon', Ulat, double(Tu(:,:,1))); shading interp; % plotting global Uwind stress
c=colorbar;
set(get(c,'title'),'string','N/m^2','fontsize',20);
%title('Global surface U-wind stress','fontsize',12);
xlabel('Longitude','fontsize',20);
ylabel('Latitude','fontsize',20);
%%Calculating Uwind stress gradient
[Tux,Tuy]=gradient(Tu); 
Tux=Tux./repmat(dx, [1, 144, 841])./(f.*P); % Incorporating the coriolis parameter and sea water density
Tuy=Tuy./dy./(f.*P);
figure
pcolor(Ulon', Ulat, double(Tuy(:,:,1))); shading interp; % plotting global Uwind stress
c=colorbar;
set(get(c,'title'),'string','N/m^2','fontsize',20);
%title('Global surface U-wind stress gradient','fontsize',12);
xlabel('Longitude','fontsize',20);
ylabel('Latitude','fontsize',20);
%%Calculating global Vwind stress
Vsq = Vwind.*abs(Vwind); % Vwind speed squared
Tv=Pa*Cd*Vsq; %Tv=global Vwind stress
Tv=permute(Tv,[1 2 3]);
figure
pcolor(Vlon', Vlat, double(Tv(:,:,1))); shading interp; % plotting global Vwind stress
c=colorbar;
set(get(c,'title'),'string','N/m^2','fontsize',20);
%title('Global surface V-wind stress','fontsize',12);
xlabel('Longitude','fontsize',20);
ylabel('Latitude','fontsize',20);
%%Calculating Vwind stress gradient
[Tvx,Tvy]=gradient(Tv);  
Tvx=Tvx./repmat(dx, [1, 144, 841])./(f.*P); %Incorporating the coriolis parameter and sea water density
Tvy=Tvy./dy./(f.*P);
figure
pcolor(Ulon', Ulat, double(Tvx(:,:,1))); shading interp; % plotting global Uwind stress
c=colorbar;
set(get(c,'title'),'string','N/m^2','fontsize',20);
%title('Global surface V-wind stress gradient','fontsize',12);
xlabel('Longitude','fontsize',20);
ylabel('Latitude','fontsize',20);
%%Calculating Ekman velocity We
Wekman = (Tvx- Tuy);
figure
pcolor(Ulon', Ulat, double(Wekman(:,:,1))); shading interp;
c=colorbar;
set(get(c,'title'),'string','m^2/s','fontsize',20);
%title('Ekman velocity for the Arctic region','fontsize',12);
xlabel('Longitude','fontsize',20);
ylabel('Latitude','fontsize',20);
ylim ([60 85]); %focusing (zooming in) on Arctic region ie. latitude above 70 degrees
end

Thanks in advance for any help, Blake.

Answer 1

[Fabrice Lambert]

lat = -90:90;
lon = -180:2:180;
Z = peaks(181);
contourf(lon,lat,Z);
hold on 
C = load('coast');
plot(C.long,C.lat,'k')

[Fabrice Lambert]

In case you don't want to plot the whole world:

xlim([min(lon) max(lon)])
ylim([min(lat) max(lat)])

will cut the coast plot to the requested range.

[Fabrice Lambert]

Update for 2021:

lat = -90:90;
lon = -180:2:180;
Z = peaks(181);
contourf(lon,lat,Z);
hold on 
C = load('coastlines');
plot(C.coastlon,C.coastlat,'k')

[Marin Vojkovic]

How would this work for just a small reagion? I have data I want to pllot over coastlines but it's just [-80 30]longitude and [30 70]latitude.

With your code I just get the whole world with a small region contourf data.

[Fabrice Lambert]

See my comment of 10 June

lat = 30:70;
lon = -80:30;
Z = peaks(111);
contourf(lon, lat, Z(1:41,:))
hold on 
C = load('coastlines');
plot(C.coastlon,C.coastlat,'k')
xlim([-80 30])
ylim([30 70])