Question about using Mapping Toolbox to project ETOPO1 data from structure array

18 views (last 30 days)
Nathan on 14 Jul 2011
Hi!
I'm currently trying to project a structure array containing 3 Fields: latitude, longitude, and altitude (bathymetry), which are 1801 x 1 double, 21601x1 double, and 1801 x 21601 int16 respectively.
The data covers latitude from -30 to 0 and longitude from 0 to 360 in 1-minute measurements (though I have a separate copy of the data with longitude measuring -180 to 180). I chose the 0-360 format because the area I need to map straddles the longitude line opposite the Prime Meridian (off the eastern coast of Australia).
I can currently plot the data in a nice-looking pcolor plot, however the map is not useful unless I can display it as a map projection (because the spacing of latitude is not uniform as the pcolor depicts but decreases as distance to the equator decreases). Ideally, I'd like to use a cylindrical mercator projection since my area of interest is near the equator.
So...after that lengthy description, how can I get my structure array containing bathymetric measurements and lat/long fields into a mercator projection? I assume the mapping tool could have some function that could help, but I'm at a loss as to how to do it. If there's a formula I could apply to the data I could do that too.
Thanks for your help, it is much appreciated!
-Nate

Mark Brandon on 31 Aug 2011
Hi Nate,
Given how long ago you asked I expect you have already got your answer. I had a similar issue with irregular lat with the Smith and Sandwell bathymetry (satbath function).
I wanted to make a large scale plot which I think that you do as well, so to get around it I just defined a new lat axis - and then re-sampled data onto that grid. My argument being at large scale the problems caused by the interpolation won't even be visible. Once on a regular grid you can just define a reference matrix and plot it anyway you like.
%%First get bathymetry
[lat,lon,elevations] = satbath(3,[-66 -48],[-71 -24]);
lon = lon - 360;
% so re-sample onto a regular grid
num_X = 350.0;
num_Y = 470.0;
change_lat = lat(:,1);
dlat = (max(change_lat) - min(change_lat))/num_Y;
lat_grid_c1 = fliplr(min(change_lat):dlat:max(change_lat))';
lat_grid_new = repmat(lat_grid_c1, 1, num_X+1);
change_lon = lon(1,:);
dlon = (max(change_lon) - min(change_lon))/num_X;
lon_grid_r1 = min(change_lon):dlon:max(change_lon);
lon_grid_new = repmat(lon_grid_r1, num_Y+1, 1);
% ZI = interp2(X,Y,Z,XI,YI)
elevations_I = interp2(lon,lat,elevations,lon_grid_new, lat_grid_new);
Then I made my reference matrix with this code
% Define a reference vector
% makerefmat Construct affine spatial-referencing matrix Syntax
%
% R = makerefmat(x11, y11, dx, dy)
% R = makerefmat(lon11, lat11, dlon, dlat)
sat_refvec = makerefmat(lon_grid_new(1,1), lat_grid_new(1,1), (lon_grid_new(1,2) ...
- lon_grid_new(1,1)), (lat_grid_new(2,1) - lat_grid_new(1,1)));
clear change_lat dlat elevations lat num_X
clear change_lon dlon lat_grid_c1 lon lon_grid_r1 num_Y
Then you can plot it as you will.
%%define bathy colourmap and plot it
bathycolourmap = [ 0.02745098 0.478431373 1 ; ...
0.176470588 0.556862745 1 ; ...
0.301960784 0.576470588 0.956862745 ; ...
0.392156863 0.670588235 0.996078431 ; ...
0.529411765 0.807843137 0.992156863 ; ...
0.529411765 0.807843137 0.992156863 ; ...
0.701960784 0.898039216 0.988235294 ; ...
0.776470588 0.941176471 0.988235294 ; ...
0.968627451 0.984313725 1 ; ...
0.980392157 0.901960784 0.756862745 ];
figure
axesm('MapProjection','mercator', 'MapLatLimit',[-64 -50], ...
'MapLonLimit',[-70 -25]);
geoshow(elevations_I,sat_refvec, 'DisplayType','surface')
colormap(bathycolourmap)
getm(gca);
setm(gca,'LabelUnits','degrees')
setm(gca,'grid','on','MLineLocation',10,'PLineLocation',10);
setm(gca,'MeridianLabel','on','ParallelLabel','on','MLabelParallel','south');
setm(gca,'MLabelLocation',10.0,'PLabelLocation',10);
set(gcf, 'Color', 'w');
tightmap
framem
And so you can see my area of interest! (I have to do something silly to plot the coastlines – which is why I am looking here).
It sounds like you may have to do a little bit with re-defining the x axis – but good luck.
Mark
Mark Brandon on 31 Aug 2011
I thought you would have moved on - you *just* want it done don't you. But glad you think the code is of use. The colourmap is one I made from the gebco scale.
Keep asking the clear Q's, and good luck, Mark