Build Lambert Grid

[carlosag]

Hello everybody,
Hope you all are fine.
I have been putting a lot of effort in learning how to construct grids to be used by ROMS.
Now, I have learnt to use both pyroms and GridBuilder but I have a question. I have noticed that the examples I have tried to use to learn these techniques are close to the equator.
I must, however, start building grids for 60˚N up to 70˚N and that would require Lambert projection. Do you recommend any particular software? I'm guessing GridBuilder might not be the best option.
Thank you in advance

[wilkin]

I have been making simple rectangular Lambert grids, optionally rotated, using simple Matlab code and the m_map toolbox and some utilities from the myroms Matlab toolset to compute the grid metrics.

Code: Select all

mstruct = defaultm('lambertstd');
mstruct.mapparallels = 71.42;
mstruct.origin = [71.42 -52.6];
mstruct = defaultm(mstruct);

% grid dimensions
Lp = 3*2^7+2;
Mp = 3*2^5+2;

% Build rectangular grid from central point
[I,J] = meshgrid(1:Lp,1:Mp);
I = I-Lp/2;
J = J-Mp/2;

% guess at suitable domain size in projected coordinates to get dx,dy 
% (alternative would be to experiment with minvtran calls to relate actual
%  distance in latitude to dy units)
XL = 0.014;
EL = 0.0035;
dx = XL/Lp;
dy = EL/Mp;

% convert i,j to x,y 
X = dx*I;
Y = dy*J;

theta = 0*pi/180;
if theta ~= 0
  % rotate the grid
  R = [cos(theta) -sin(theta); sin(theta) cos(theta)];
  XY = [X(:) Y(:)]';
  XY2 = R*XY;
  X = reshape(XY2(1,:),size(X));
  Y = reshape(XY2(2,:),size(Y));
end

% convert projected coordinate back to lon/lat
[rlat,rlon] = minvtran(mstruct,X,Y);

% Coastline for plotting
% clear S
if ~exist('S','var')
  file = fullfile('/Users/wilkin/Dropbox/MATLAB/toolbox/myroms/m_map',...
    'private/gshhs_h.b');
  indexfile = gshhs(file, 'createindex');
  S = gshhs(file,range(rlat),range(rlon));
end

% Plot

figure(1)
clf
axesm(mstruct,'MapLonLimit',range(rlon),'MapLatLimit',range(rlat),'Grid','on')
plotm(rlat,rlon,'k.')
a = axis;
levels = [S.Level];
L1 = S(levels==1);
geoshow([L1.Lat],[L1.Lon],'Color','k')
axis(a)

figure(2)
clf
plot(rlon,rlat,'.')
a = axis;
hold on
for k=1:size(S)
  plot(S(k).Lon,S(k).Lat,'-')
end
hold off
axis(a)
amerc

% write the roms grid file
set_grid(rlon, rlat, 'test.nc','linear','i')
g = roms_get_grid('ks2.nc');

[carlosag]

Dear Wilkin,
Thank you very much.
I will give it a try.
Does the m_map tool also helps to smoothing the bathymetry and helps with the masks as GridBuilder does?

[wilkin]

I misspoke in my message. I previously used the contributed m_map mapping toolbox to do this, but the code I posted uses Matlab's own Mapping toolbox ... the mstruct = defaultm('lambertstd'); etc steps.

It's just very simple code to make a simple grid. In this example we created the bathymetry and land/sea mask separately, which I generally prefer to do to have control over smoothing and paying attention to sill depths, etc. That said, GridBuilder usually does a pretty good job for grids with resolutions that its database is suited to.

In my example code I posted, this is a small fjord in West Greenland and the global bathymetry databases are definitely not adequate.

[carlosag]

Thank you very much, I will try your suggestion and try to understand it!!
All best