ROMS.jl/src /create_grid.jl produces incorrect x_rho, y_rho, x_psi, x_u, x_v, etc.

[lnferris]

the spacing of grids x_rho, y_rho, x_psi, x_u, x_v is inconsistent (e.g ~20 times smaller) than that of lon_rho, lat_rho, etc.

[lnferris]

This is not very elegant but I believe this is a pseudocode (matlab) solution that preserves the grid spacing of lon_rho, lat_rho, etc. Please check that you agree:

[x_rho,y_rho] = map_to_grid(lon_rho,lat_rho,0,0);
[x_u,y_u] = map_to_grid(lon_u,lat_u,0.5,0);
[x_v,y_v] = map_to_grid(lon_v,lat_v,0,0.5);
[x_psi,y_psi] = map_to_grid(lon_psi,lat_psi,0.5,0.5);

function [x,y] = map_to_grid(lon,lat,xshift,yshift)

    dlat_m = 110.574e3;
    dlon_m = 111.320e3*cosd(lat);
    
    y = zeros.*lat;
    for j = 2:size(lon,2)
     y(:,j) = dlat_m.*(lat(:,j)-lat(:,j-1))+ y(:,j-1);
    end
    y = y + yshift.*(y(:,2)-y(:,1));
    
    x = zeros.*lon;
    for i = 2:size(lon,1)
        x(i,:) = dlon_m(i,:).*(lon(i,:)-lon(i-1,:))+ x(i-1,:);
    end
    x = x + xshift.*(x(2,:)-x(1,:));

end

[Alexander-Barth]

I just fixed the scaling issue. Reading thought the seagrid code, it seems that these x/y are expressed in the Mercator projection (or other conformal projection).

Apparently, in the "official" c_grid.m, these variables x_/y_ are no longer written for spherical grids
https://github.com/myroms/roms_matlab/blob/main/grid/c_grid.m

I just tested that indeed ROMS runs just fine without these variables.

In future we might decide to just drop these variables for spherical grids.

[lnferris]

Okay. I think the scaling is still incorrect.

[Alexander-Barth]

Can you explain what you mean?
There is no map projection without distortion. At least the limitations of the Mercator projection are well understood.

Your proposed code gives the following grid:

For the code:

lonv = 0.:2:50
latv = 0.:2:60

lon = [x_ for x_ in lonv, y_ in latv]
lat = [y_ for x_ in lonv, y_ in latv]


xshift = yshift = 0

dlat_m = 110.574e3
dlon_m = 111.320e3 .* cosd.(lat)

y = 0 .* lat;
for j = 2:size(lon,2)
   y[:,j] = dlat_m .* (lat[:,j]-lat[:,j-1]) +  y[:,j-1];
end

y = y .+ yshift .* (y[:,2]-y[:,1]);

x = 0 .* lon;
for i = 2:size(lon,1)
   x[i,:] = dlon_m[i,:].* (lon[i,:]-lon[i-1,:]) + x[i-1,:];
end

x = x .+ xshift * (x[2,:]-x[1,:])'

using PyPlot
clf()
plot(x,y,"k-")
plot(x',y',"k-")

The projected grid is no-longer orthogonal. Even if the distance between direct cell neighbors (along x and y axis) is correct in your projection, the distance between e.g. the points i,j and i+1,j+1 is not.

The Mercator projection preserve angles and but distances have to be interpreted using a scaling function.

For reference here is what pyroms does:

https://github.com/ESMG/pyroms/blob/1370d22ab94c877648f82ffc440d25d782dfce62/pyroms/pyroms/sta_hgrid.py#L57-L85

The projection can be chosen by the user. I think it will implemented that.

(it would be easier for me if you can just share julia code in future, thanks!).

[lnferris]

I was comparing dx and dy between adjacent points to 1/pm and 1/pn to check equivalence. I did not realize most users optimize for planar orthogonality and apply a scaling function to the data to recover the correct scale.

I tried to implement a preservation of scale using these two modifications:
create_grid.jl
projections.jl

which resulted in this grid that preserves the distance found between points on the spherical grid. I think it is okay due to spherical orthogonality (e.g. #11) but could be wrong.