Surface geostrophic velocity issue

[donoso]

Hello guys,

I'm trying to do some post processing of a curvilinear application, such as computing the surface geostrophic velocity. To do that, I'm derivating zeta field along it's cols and rows (which correspond to xi and eta coordinates), dividing by it's metric resolution (1/pm, 1/pn) and then multiplying bythe term g/f. Finally, I transformed the vector components from xi and eta to lon and lat coordinates by means of:

ug[lon,lat] = ug[xi,eta].*cos(angle) - vg[xi,eta].*sin(angle);
vg[lon,lat] = vg[xi,eta].*cos(angle) + ug[xi,eta].*sin(angle);

However, as you may see in the attached figure, the resulting vector components [lon,lat] are not tangential to the zeta isolines at the lower right corner of domain. So, do you know what could be the issue?

Thanks in advance!

Cheers,
David

[wilkin]

pm, pn have units inverse distance, so d(zeta)/d(xi) is diff(zeta) TIMES pm.

But pm,pn are defined at cell center rho points, same as zeta. What ROMS does is average pm to the cell face u-point, say pm_u, and calculate diff(zeta)*pm_u. Then you have the zeta slope at u-points, and you need to average them back to the zeta (rho) point to apply the angle rotation. Likewise for d(zeta)/d(eta) and pn an v component.

At the perimeter of the domain, the step of the average back to rho is going to leave you a row and column short, so you need to pad these by duplicating the lost row/columns.

[donoso]

Thank you John by your feedback!

I think we may do the overall calculus at different order, but the results will not be very different. By the way, I made de calculus as you said according the following MATLAB code lines:

% Spatial derivates (adimensional)
dzetadxi = diff(zeta,1,2) ./ rho2u_2d(1./pm); % at U-points
dzetadeta = diff(zeta,1,1) ./ rho2v_2d(1./pn); % at V-points

% Geostrohic velocities (m/s) >> [Xi,Eta] coordinates
vgeta = +(g./f) .* u2rho_2d(dzetadxi); % at RHO-points
ugxi = -(g./f) .* v2rho_2d(dzetadeta); % at RHO-points

% Rotate vector field from [Xi,Eta] to [Lon,Lat]
[ug,vg] = rotate_vec(ugxi,vgeta,angle,1); % at RHO-points

Where;
- "rho2u_2d.m/rho2v_2d.m" are Pierrick Penven's (ROMS-Agrif) functions which transform a field at rho points to a field at u/v points).
- "u2rho_2d.m/v2rho_2d.m" are Pierrick Penven's (ROMS-Agrif) functions which transform a field at u/v points to a field at rho points).
- "rotate_vec.m" is a Hernan G. Arango's function (ROMS-Rutgers) which rotate vector field from [Xi,Eta] to [Lon,Lat]).

I attached here a new plot with the results (considering an extended domain at eta direction to better representation of distances ratio). As you may see, the issue persists. I refuse to think that something is wrong with the angle variable, so I'm thinking about eventual problem with the visualization... I'm suspicious of "quiver.m" function (which do vectors representation).

Cheers,
David.

[wilkin]

I think your average to u-points of 1/pm is not the same as 1/(average of pm).

ROMS uses (look in metrics.F)

Code: Select all

!-----------------------------------------------------------------------
!  Compute m/n, 1/m, and 1/n at horizontal U-points.
!-----------------------------------------------------------------------
!
      DO j=JstrT,JendT
        DO i=IstrP,IendT
          pmon_u(i,j)=(pm(i-1,j)+pm(i,j))/(pn(i-1,j)+pn(i,j))
          pnom_u(i,j)=(pn(i-1,j)+pn(i,j))/(pm(i-1,j)+pm(i,j))
          om_u(i,j)=2.0_r8/(pm(i-1,j)+pm(i,j))
          on_u(i,j)=2.0_r8/(pn(i-1,j)+pn(i,j))
        END DO
      END DO

[donoso]

I'm agree the expression "rho2u_2d(1./pm)" is not exactly the same as "1./rho2u_2d(pm)" [which, according to your comments, seems to me to be equal to "om_u(i,j)=2.0_r8/(pm(i-1,j)+pm(i,j))" at metrics.F]. But that difference is quite small and could not explain the presented issue.

I attached here the content of rho2u_2d.m function:

Code: Select all

function var_u=rho2u_2d(var_rho);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% function var_u=rho2u_2d(var_rho);
%
% interpole a field at rho points to a field at u points
%
% input:
%
%  var_rho variable at rho-points (2D matrix)
%
% output:
%
%  var_u   variable at u-points (2D matrix)  
% 
%  Further Information:  
%  http://www.croco-ocean.org
%  
%  This file is part of CROCOTOOLS
%
%  CROCOTOOLS is free software; you can redistribute it and/or modify
%  it under the terms of the GNU General Public License as published
%  by the Free Software Foundation; either version 2 of the License,
%  or (at your option) any later version.
%
%  CROCOTOOLS is distributed in the hope that it will be useful, but
%  WITHOUT ANY WARRANTY; without even the implied warranty of
%  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%  GNU General Public License for more details.
%
%  You should have received a copy of the GNU General Public License
%  along with this program; if not, write to the Free Software
%  Foundation, Inc., 59 Temple Place, Suite 330, Boston,
%  MA  02111-1307  USA
%
%  Copyright (c) 2001-2006 by Pierrick Penven 
%  e-mail:Pierrick.Penven@ird.fr  
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[Mp,Lp]=size(var_rho);
L=Lp-1;
var_u=0.5*(var_rho(:,1:L)+var_rho(:,2:Lp));
return

[wilkin]

I added a function roms_surface_geostrophic_velocity.m to my ROMS tools at https://github.com/johnwilkin/roms_wilkin. It returns surface geostrophic velocity in east-north coordinates on the lon_rho grid ready for plotting by quiver. It gives vectors that align with zeta contours in my test plots.

This works in the same framework of my other tools, so you'll have to read the grd structure with roms_get_grid.m etc., but I think it is sufficiently well documented internally you'll get it to work. Let me know if it gives the same or different results to your own code. I recycled a lot of code chunks that are in other mature functions in that repo, so I'm pretty confident it's correct.

[donoso]

Thank you John, I appreciate your effort. I downloaded your code but the results are almost the same, the only differences are velocity magnitudes due to used precision at gravity of Earth.

On the other hand, the complete domain of my application include West Antarctic Peninsula. In this sense, I attached a new plot which shows a bigger area and also a greater range of the angle. As you may see, at [xi,eta] coordinates (a and b panels) the velocity components (subsampled each 5 grid points) are aligned with zeta contours. But at [lon,lat] coordinates (c and d panels) the rotated velocity components are not tangential to the zeta isolines both southward and northward to South Shetland Islands among others locations. In this moment I'm a little perplexed, so do you think the issue would be atribuible to something wrong with the angle array?

Cheers,
David.

[wilkin]

Yes, I would suspect your grid angle.

How did you generate the grid? Was it an angle preserving map projection?

Try plotting some velocity vectors with my roms_quivergrd and look closely if they are parallel to the coast. They absolutely should be.

[msd]

Hi John,

This grid is based on an earlier grid made at a coarser resolution using the old Fortran gridpak package (so the grid angle would have come from "sphere"). I didn't make the finer resolution grid, so I don't know how the grid angle was changed, but I'm guessing some simple interpolation that shouldn't have changed things much.

Note that the angle of error of the vector w.r.t. the isolines of zeta seem to be more a function of the orientation of the isolines than where in lat/lon space the vector is (e.g. vectors that are due N or S seem to be correct no matter what longitude, and maybe vectors that are due E or W are OK too). My guess is this shows that it may not be an issue with the model "angle" variable since that varies noticeably by longitude over this bigger plot.

Not that I have a good other idea for what's going on...

[arango]

Use Matlab script grid/roms_metrics.m to recompute the angle and other quantities to see what you get. You can download that script from the ROMS Matlab repository. Then, compare your rotation angle arrays. Interpolation is not a good idea!

[wilkin]

Well, this is a mystery.

I recomputed your angle and metrics from the grid lon/lat using roms_metrics.m (part of the Matlab tools at myroms.org). The values are not appreciably different.

So, I made a sanity check plot using mature code that that I have used for years. I plotted two sets of vectors with u=0 and v=1, and with u=1 and v=0. These should plot in line with the grid coordinates. They do not. Here, with u=0 and v=1 in green, and u=1 and v=0 in red. (By this I mean the u,v were placed on their native Arakawa-C grid, and my roms_quivergrd code averaged them to the rho point cell centers and applied the angle rotation).

Could it be that at very high latitude something is awry with the angle calculation? If so, that potentially means that all results with this grid have misaligned the direction of winds if you use the automatic regridding option.

[msd]

Luckily, I pre-interpolate and rotate the winds and do not use the automatic regridding option in ROMS.

Still strange though: I don't use MATLAB to make vector plots, so I haven't used quivergrd (or quiver) at any latitude. I have made plenty of vector plots of model velocity fields at high latitudes where the grid (and vector fields) were rotated back to true N. I've definitely had examples where model velocities followed along bathymetry quite nicely (at least to the detail I was looking at) as expected along steep slopes where there wasn't anything to break potential vorticity constraints.

[arango]

Hmm, the grid looks very orthogonal to me. It is hard to think that it is due to the averaging of the quiver plot. The angle between the red and green vectors is pretty much constant and around 30 degrees. It is weird. I wonder what do we get with uv_rotate.F when we output the vectors to geographical East and Northward directions.

[msd]

Yeah, I don't get why the angle between the red and green vectors isn't 90 degrees if the grid looks pretty orthogonal and the centers are at close to right angles on the screen (at least it looks that way on my screen). If the ROMS "angle" was off, I could see why the arrows might not line up with the grid. However, if you're doing a simple rotation of [1,0] and [0,1], I would expect the two vectors to still be normal to each other unless there was some issue with the plotting of the vectors on the page? Or am I missing something?

[donoso]

I attached here a new plot that shows an example of model output's bottom velocities (subsampled each 2 grid points), including bathymetry contours every 100 m up to 500 m depth (dashed lines) and every 500 m up to 5000 m depth (solid lines). As Mike said, I would expect that velocities would follow the bathymetry over steep slope but that is not very clear to me over Antarctic Slope and northward to Bransfield Strait.

[arango]

Are you expecting flow to be everywhere along bathymetric contours? We can have bottom boundary layers, turbulence, pressure gradient differences, and friction that may force the flow across isobars and bathymetric contours. I think that what you need to plot are the barotropic currents at every grid point, and check how that looks. Now, when drawing vectors, you need to consider the vector's scale and what part of the vector is located at the (lon, lat) grid cell. The tail, center, or head of the vector? It seems to me that making conclusions on a vector plot is not that robust.

[donoso]

OK Hernan, I made a new plot with barotropic velocities at every grid point. The vectors are plotted using "quiver.m" at the same fashion of "roms_quivergrd.m", I included a vector's scale over the Antarctic Peninsula and I think vector 's tails are located at [lon,lat] grid cells.

[msd]

I really like John's sanity check, so I did the same thing with what I use to plot vectors.

Each vector pair in the image is either [1,0] or [0,1] along the model xi- and eta- axes, rotated using the model "angle" variable, and the vector tails start at the rho-point. This is using every fourth model grid point. The code that draws the vectors is pretty simplistic and needs explicit start and end grid points (in lat/lon space in this case). Therefore, I have to figure out a ratio between distance along the east axis and distance along the north axis so that the vector angles get drawn correctly. As a sanity check of that, the bottom right has four vectors all with tails starting in the same location and values along the East/North axes: [2,0], [0,2], [sqrt(2),sqrt(2)], and [2,2].

Note that the physical scale along the x and y-axes (km along the earth/inch on the page) is not the same (although it's close with these boundaries, especially at the bottom of the plot), but the scaling that I use between east/north is set so that a vector of a given magnitude will have a constant length on the page no matter what angle it's pointed at.

Everything looks pretty reasonable to me in this plot. I could be messing up how I do the scaling between east/north and just getting lucky (i.e. making the same mistake that the angle generating code is making). However, could there be some issue with how quiver draws vectors at high latitudes? Again, I don't use MATLAB to make vector plots, so I don't have any experience with this.

This was bugging me (and it's too hot outside here to do the yard work I should be doing), hence why I'm making pictures on a Saturday...

[wilkin]

Mike is on to the problem. It's a failing of Matlab quiver with no easy workaround.

If the plot space has a highly elongated aspect ratio the quiver vector components are distorted. This is why one can't simply put vectors of horizontal and vertical velocity on a vertical slice plot.

I have a rather clunky function called roms_curquivergrd.m

Code: Select all

function [han,data] = roms_curquivergrd(u,v,grd,lon0,lat0,td,nd,...)

... that draws what are essentially particle trajectories of duration td days initialized at a set of requested starting locations (lon0,lat0). I repeated my unit vector test with roms_curquivergrd in the attached plot. The red and blue particle trajectories for x and y-dir unit vectors now align with the coordinate grid as they should.

roms_curquivergrd.m works by normalizing the velocity components by the grid metric factors, pm and pn, to get velocity in fractional i,j per second, then multiplying by the requested duration (td) to get final i,j locations, then maps the result back to lon/lat space for plotting.

It's a bit slow because it does a 4th order Runge-Kutta calculation for the trajectories, which is probably overkill, so choosing a large number of lon0,lat0 starting positions can be very slow. But the lon0,lat0 are arbitrary 1-D vector so can be chosen carefully to illustrate the flow field without excessive compute burden. The code is part of https://github.com/johnwilkin/roms_wilkin

[arango]

Great detective work! I also was starting to suspect about the quiver function or its usage. I am planning to draw the vectors with ROMS native NCAR graphics software next week to see what I get. I am wondering what kind of results we will the m_map (m_quiver) tool in Matlab will give. The vectors are supposed to be scaled, but maybe we get the same result as the default quiver function.

[johnluick]

A very interesting thread. Am I correct in thinking that it implies Matlab's map-viewed quivers are in error by about 15% at around 60 degrees latitude (and getting worse as latitude increases)?

I happened to have this article on my screen at the same time while reading down through it - it purports to explain why it was too hot for David to do yard work. https://edition.cnn.com/2021/08/26/worl ... index.html

[donoso]

Thank you for your comments and work!

I was doing some tests with vector plots over a map projection and I think could be useful to share it. I attach two Figures with the first areas that I showed (at the beginning of this post) be part of whole West Antarctic Peninsula domain. The left panels shows the geostrophic velocities at (unprojected) [lon,lat] coordinates while the right panels shows the geostrophic velocities at projected Cartesian coordinates over Mercator ([lon,lat]*). As we can see, the projected velocities (right panels) are now aligned with zeta contours, which is great to me and probably Mike's too. Although it is nice to know that this problem is only a visualization matter, I'm still wonder if it could be present at curvilinear applications only?

Cheers,
David.

NOTE: I made all vector plots with quiver.m but on the projected cases additionally I used mfwdtran.m (from Mapping Toolbox -a native MATLAB library-), which transforms unprojected [lon,lat] coordinates to the projected Cartesian coordinates frame using the map projection defined by the user.

[arango]

Great, I think that we need to design a Matlab function that plots vectors correctly regardless if the grid is curvilinear or not.

[rduran]

When everything else fails, read the manual. From Matlab's quiver documentation:

By default, the quiver function shortens the arrows so they do not overlap. Call axis equal to use equal data unit lengths along each axis. This makes the arrows point in the correct direction.

A simple

Code: Select all

axis equal

should do it.

[donoso]

Hi Rodrigo,

I think the main problem is not the proportionality (aspect ratio) between the coordinate axes. If that were the case, the appearance of the zeta contours would also be affected so that they could never be aligned with the velocity vectors.

[rduran]

I see your point, but reading he manual still applies since reading the help for m_quiver I am directed to m_vec for this type of problem. Hope that helps.

[wilkin]

David is correct. Setting

Code: Select all

>> axis equal

does not solve the problem. Try it and you'll see.

[rduran]

I know, that is why I mentioned that I see his point, and that m_vec may solve the problem. in m_quiver help it says

Note that we do not scale arrows with respect to map
coordinates! Instead, the arrows will correspond better to actual motions
over some time step. The tradeoff is that a single scale arrow cannot
be accurate for the entire map (M_VEC scales arrows according to
map coordinates).