Task: Extract a high resolution elevation data (etopo2) covering Chesapeake Bay. Create a map and contour the elevation/depth for this region.
MATLAB contains a variety of files containing geographic and political information as well as DEM (digital elevation maps). There is too much information to cover in detail. The files below provide (mainly) elevation information although one of them provides coastline values. We have already used 'landareas.shp' providing land areas as a GIS shape file.
The built-in topo data, available with load topo provides 1 degree by 1 degree (180x360) elevation and depth information. The elevation is in variable topo. Along with it comes two colormaps (topomap1 and topomap2) as well as latitude and longitude range information.
A coarse resolution coastline data base is available with load coast which produces variables lat and long.
A more convenient way to use this data is to read it to a structure with coast=load('coast'). The variables are now coast.lat and coast.long, which means that they are not likely to conflict with other active variables.
The above elevation data set (topo) is a bit crude, with a value every degree of longitude and latitude. A finer data set is available with tbase which reads a DEM (etopo5) with a resolution of 5 minutes (or 1/12th of a degree). The command reads the tbase.bin data file (use the link to get your copy of this data file. Note: it is 18Mb). The following command will extract a 1 degree elevation data set from etopo5 for the indicated latitude and longitude ranges.
SF=12; % subsample factor (every twelveth value) LatLim=[0 20];LonLim=[-100 -60]; [Elev, RVec]=tbase(SF,LatLim,LonLim)Make sure that the data file is in a folder in your current path so MATLAB can find it.
An even finer elevation data set is available at 2 minutes (1/30 th of a degree) with satbath. It requires the topo_8.2.img data file (use the link to get your copy. It is 131 Mb). A similar extraction command is
SF=10; % subsample factor (every tenth value) LatLim=[0 20];LonLim=[-100 -60]; [LAT LON,Elev]=satbath(SF,LatLim,LonLim);A contour plot of this data is produced with contour(Lon,Lat,Elev) or it can be plotted on a map projection with contourm.
There is a very high precision DEM (gtopo30) which is a 30 second (half minute or 1/120 th of a degree) DEM. This data is stored in sectors which makes extraction of the data a bit complicated. The routine to extract this data is gtopo30. Read the help information for details on the use of this data. You will also need to download the data file; information on doing this is in the MATLAB documentation.
Flow chart for task
%%% Set up location for chesapeake bay %%% extract DEM for area %%% get size of extracted DEM %%% create Lat and Lon arrays for grid locations %%% create base map %%% contour the DEM %%% add lat,lon tick marks
Script for task
%%% Set up location for chesapeake bay
SF=1; % subsample factor (every value)
LatLimCB=[36.5 39.5];LonLimCB=[-77 -75];
%%% extract DEM for area
[Elev, RVec]=tbase(SF,LatLimCB,LonLimCB);
%%% get size of extracted DEM
[nr nc]=size(Elev);
%%% create Lat and Lon arrays for grid locations
LatCB=linspace(LatLimCB(1),LatLimCB(2),nr);
LonCB=linspace(LonLimCB(1),LonLimCB(2),nc);
%%% create base map
figure
axesm('MapProjection','eqaconic',...
'MapLatLimit',LatLimCB,'MapLonLimit',LonLimCB,...
'Grid','on','Frame','on')
%%% contour the DEM
[con han]=contourm(LatCB,LonCB,Elev,[-80:5:10]);
% clabelm(con,han)
%%% add lat,lon tick marks
mlabel ON;plabel ON;