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

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

Answer 1

[Mark Brandon]

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.

So what I did for my plot was this.

%%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

[Nathan]

This is what we ultimately decided to do, though I ended up switching to GMT to do the mapping side of things. Thanks for the detailed code breakdown, I do still need to do some stuff with MATLAB, so this will really help me a TON in a week or two!

[Mark Brandon]

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