EASYGRID: Grid generation for ROMS

[Diego]

EASYGRID is a simple matlab script to make ROMS grids. It is an austere alternative, but it fairly easy to get it up-and-running (needs MEXNC, SNCTOOLS and ROMS-Matlab-Toolkit). This Matlab script also allows the user to (1) smooth bathymetry, (2) interactively edit land/sea mask, and (3) create simple initialization NetCDF files for ROMS.

HERE is a tutorial that explains how to download, install and use EASYGRID.

I see EASYGRID as a beginner's tool, yet capable of making grids for realistic applications. Its simplicity may make it appropriate for learning (and teaching) the steps required to create a grid for ROMS. It should also facilitate new users to start experimenting with their own projects in ROMS... hopefully speeding up the learning process.

Comments, questions, bug-reports are all welcome in this thread.

I hope this helps.

Diego

[lindseyk]

Thanks for EASYGRID, it's very helpful to see the steps that go into grid generation in explained and executed in one script. The tutorial is also extremely helpful, especially including step by step instructions for including mexnc and snctools, which I was never quite sure I'd done right. I'm still a novice at grid development but have recently managed to produce a grid using SeaGrid, but wanted to sea how this would compare and to get tips on writing gridded initial condition files.

So I decided to try out easygrid to generate a very large southern ocean grid. Althought it takes a while to run, things go fine until I get the following error:

GENERATING grid...
??? Error using ==> qhullmx

While executing: | qhull d Qt Qbb Qc

Options selected for Qhull 2003.1 2003/12/30:
delaunay Qtriangulate Qbbound-last Qcoplanar-keep _pre-merge
_zero-centrum Pgood Qinterior-keep _max-width 3.6e+02
Error-roundoff 5e-13 _one-merge 3.5e-12 Visible-distance 1e-12
U-coplanar-distance 1e-12 Width-outside 2e-12 _wide-facet 6e-12

Last point added to hull was p903479.
Last merge was #3497178.

Qhull has finished constructing the hull.



At error exit:



Error in ==> delaunayn at 114
t = qhullmx(x', 'd ', opt);

Error in ==> griddata>linear at 151
tri = delaunayn([x y]);

Error in ==> griddata at 120
zi = linear(x,y,z,xi,yi,opt);

Error in ==> easygrid at 323
h = griddata(xbathy,ybathy,zbathy,lon_rho,lat_rho,'linear');


Does his have to do with the size of my grid? These are it's dimensions, I'd like to work more on the kilometer scale:

X = 7900000; % Width of domain (meters)
Y = 8300000; % Length of domain (meters)
rotangle = 0; % Angle (degrees) to rotate the grid
resol = 50000; % Cell width and height (i.e. Resolution)in meters.
N = 23; % Number of vertical levels

Thanks!

[Diego]

Hi lindseyk,

Sorry for the late reply... I am just coming back from a conference trip.

I tried Easygrid with your dimensions (listed above) and it worked, but the grid was all funky... Easygrid cannot handle grids over SUPER LARGE domains. However, if I take one zero off X, Y and resol (i.e. make the grid 10 time smaller)... Easygrid works fine.

Cheers!

[turuncu]

Hi,

What about lakes? I try to generate grid for lake and it gives following results (see attached image). I also faced some problems defining X,Y and resol parameter. I know the size of the lake (122 km, 81 km) and i set the resol parameter as 1000 m. So in this case i expect that grid must be cover whole lake but it is not. I plot the image using parameter as X = 700000; Y = 500000; and resol = 10000;.

Best,

--ufuk

grid plot

Picture 1.png (64.38 KiB) Viewed 45068 times

[Diego]

Hi ufuk,

I am not sure I understand your problem...

For what I understood, you estimated X Y and resol and when you plug your estimated values in Easygrid, the resulting grid-plot did not look as you expected... it was too small or not quite covering the whole lake. Did I get it right?

If so... you may have to iteratively change (one variable at the time) X, Y, lat, lon, resol and/or rotangle... plot in between to see the effect of each change you made. I know it is a rather primitive way of adjusting a grid, but in ~10 minutes you should be finished. Not elegant, but it works.

Diego

[eastwind111]

When I try to use fillseg(coast) to get a coastline, some error appeared:
???Input argument "c1" is undefined.

Error in ==>fillseg at 20
if(length(c1(:))==3)

The usage of this function is [h]=fillseg(x,c1,c2),
but it seems no problem for some others when input "(coast) "

[ocecept]

Sorry Diego, In my previous message, I was supposed to answer the question not repeat it.

I had the same problem reported by user "eastwind11" and the solution as Diego said is:

fillseg(coast_new,[0.5 0.5 0.5],[0 0 0])

It will result gray in land areas and black coastline.

First, make sure you have jointed the coast segments using:

coast_new=join_cst(coast,.0001);


Cheers;


_________________
Carlos Teixeira
PhD Student
SARDI/ UNSW
SA - Australia

[Diego]

eastwind111, ocecept:

fillseg is explained in the WikiROMS page https://www.myroms.org/wiki/index.php/seagrid

...there is also a very nice movie-tutorial ("Grid Generation with SeaGRID" at: http://www.myroms.org/index.php?page=tutorials ), where the use of fillseg is explained.

The way fillseg is used in these tutorials is...

Code: Select all

>> fillseg(x);

I have only used it like that (i.e. without h, c1 and c2)...

It seems that c1 and c2 are used to define the color of land and coastline... in which case each needs to be a 3-number array (between 0 and 1) defining a RGB color... [so for example c1=[0.5 0.5 0.5]) ] ...but as I said, I haven't used fillseg with is functionality.

Cheers,

Diego

[ocecept]

Hi all.

First thanks for EASYGRID. I was able to create the grid to my area, smooth the bathymetry and mask the coast with no problems. My grid has 955 Km Width, 755 Km Length, 6 Km resolution and is rotated -33 degrees. Now, I'm using the EASYGRID default "spherical" equals "F".

I'm forcing a simple 2D model with a periodic wind and it's working fine with the expected results.

The problem is: My ROMs output files (ocean_his.nc, ocean_avg.nc, etc) don't have any of the latitude and longitude variables (lon_psi, lon_rho,lon_u,lon_v, etc).

From the wrt_info.f90 at Build folder we have:

Grid coordinates of RHO-points.
!
IF (spherical) THEN
scale=1.0_r8
wrt_info=nf90_inq_varid(ncid, 'lon_rho', varid)
IF (ncid.eq.ncSTAid(ng)) THEN
ELSE
wrt_info=nf_fwrite2d(ng, model, ncid, varid, 0, r2dvar, &
& LBi, UBi, LBj, UBj, scale, &
& GRID(ng) % rmask(LBi,LBj), &
& GRID(ng) % lonr(LBi,LBj))
END IF
IF (wrtThread.and.(wrt_info.ne.nf90_noerr)) THEN
WRITE (stdout,10) 'lon_rho', TRIM(ncname)
exit_flag=3
ioerror=wrt_info
RETURN
END IF
wrt_info=nf90_inq_varid(ncid, 'lat_rho', varid)
IF (ncid.eq.ncSTAid(ng)) THEN
ELSE
wrt_info=nf_fwrite2d(ng, model, ncid, varid, 0, r2dvar, &
& LBi, UBi, LBj, UBj, scale, &
& GRID(ng) % rmask(LBi,LBj), &
& GRID(ng) % latr(LBi,LBj))
END IF
IF (wrtThread.and.(wrt_info.ne.nf90_noerr)) THEN
WRITE (stdout,10) 'lat_rho', TRIM(ncname)
exit_flag=3
ioerror=wrt_info
RETURN
END IF
ELSE
scale=1.0_r8
wrt_info=nf90_inq_varid(ncid, 'x_rho', varid)
IF (ncid.eq.ncSTAid(ng)) THEN
ELSE
wrt_info=nf_fwrite2d(ng, model, ncid, varid, 0, r2dvar, &
& LBi, UBi, LBj, UBj, scale, &
& GRID(ng) % rmask(LBi,LBj), &
& GRID(ng) % xr(LBi,LBj))
END IF
IF (wrtThread.and.(wrt_info.ne.nf90_noerr)) THEN
WRITE (stdout,10) 'x_rho', TRIM(ncname)
exit_flag=3
ioerror=wrt_info
RETURN
END IF
wrt_info=nf90_inq_varid(ncid, 'y_rho', varid)
IF (ncid.eq.ncSTAid(ng)) THEN
ELSE
wrt_info=nf_fwrite2d(ng, model, ncid, varid, 0, r2dvar, &
& LBi, UBi, LBj, UBj, scale, &
& GRID(ng) % rmask(LBi,LBj), &
& GRID(ng) % yr(LBi,LBj))
END IF
IF (wrtThread.and.(wrt_info.ne.nf90_noerr)) THEN
WRITE (stdout,10) 'y_rho', TRIM(ncname)
exit_flag=3
RETURN
END IF
END IF
!

Witch I understand that in order to have the lat and lon variables I need to use spherical coordinates. Is that true? Do I have any implications using spherical coordinates?

I also noticed that EASYGRID don't save the ANGLE information to the grid file and I couldn't find where ROMs use it. Do you know what is the implications of don't have this information?

Many thanks;

[kate]

The idea behind the spherical flag is that you either have a grid on the sphere that has lat,lon or you have a grid in some cartesian x,y space. If it is x,y, you don't want have to invent some meaningless lat,lon values for it.

As for the angle, it is needed for instance if you ask ROMS to interpolate coarse NCEP winds onto your ROMS grid. If you do that, you need to rotate the winds to align with your grid instead of the NCEP east-west and north-south.

[Diego]

... the script (at least the version I have) does not work out the angle of the grid and write it on the nc file.

Is there any updated version of easygrid where the angle is an output in the netcdf file?

You can add the a small piece of code to easygrid so that it writes the angle to the nc file... here is how:

1) Find the next piece of code in easygrid (last entry in the SAVE GRID FILE section)

%--------------------------------------------------------------------------
dims = { 'one'};
nc{ 'dy'} = ncdouble(dims);
nc{ 'dy'}(:) = resol ./ (60*1852); % Estimated resolution in y (degrees)
nc{ 'dy'}.description = 'resolution in y (degrees)';

2) Copy-paste the following code below it :

%--------------------------------------------------------------------------
dims = { 'one'};
nc{ 'rotangle'} = ncdouble(dims);
nc{ 'rotangle'}(:) = rotangle; % angle of conterclock-wise grid rotation (degrees)
nc{ 'rotangle'}.description = 'angle of conterclock-wise grid rotation (degrees)';

If you require a specific angle variable name in the nc file (which I don't know if you do)... you can change the 'rotangle' name (that is nc{'rotangle'} ) to be whatever you want.

Check the other entries in the "SAVE GRID FILE section"... It is relatively easy to figure out the logic used to save nc varibles. That way you can add ANY matlab variable to your nc files (grid or initialization).

Hope this helps.


Diego

[wgjlove]

Hello,everyone.
I need some help~~
I use easygrid to make grid of East China Sea,117E-131E,24N-41N
First of all,I should get my own bath.mat and coast.mat
so I use ROMS matlab-toolkit(http://www.myroms.org/index.php?page=Processing) to do such things
1.for the bath\extract.m:
Llon=117.0; % Left corner longitude
Rlon=131.0; % Right corner longitude
Blat=24.0; % Bottom corner latitude
Tlat=41.0; % Top corner latitude
run and get bath.mat
2.for the coastlines\extract.m:
Llon=117.0; % Left corner longitude % East of China shelf Sea
Rlon=131.0; % Right corner longitude
Blat=24.0; % Bottom corner latitude
Tlat=41.0; % Top corner latitude
run and get coast.mat
3.for the easygrid.m:
lat = 24; % Latitude (degrees) of the bottom-left corner of the grid.
lon = 117; % Longitude (degrees) of the bottom-left corner of the grid.
X = 1420000; % Width of domain (meters)
Y = 1900000; % Length of domain (meters)
rotangle = -1; % Angle (degrees) to rotate the grid counterclock-wise
resol = 25000; % Cell width and height (i.e. Resolution)in meters. Grid cells are forced to be (almost) square.
N = 16; % Number of vertical levels
run ,but I got a terrible figure,that water didn't cover the sea but cover the land
I really appreciate your help!
perhaps u ask me why don't follow the way that "https://www.myroms.org/wiki/index.php/seagrid" tells us ,because I can't understand the following words :'This will also sort the coastline segments so that the longest segment is first. Often this segment is the main coastline, while the rest of the segments are islands. To see what has happened, clear the figure and try "fillseg" again'

[perezro]

Hello!

I just sort the problem with the coastline that I believe Is the error you are in the middle of...
After getting the coastline most of the times I needed to re-adapt a bit the data so you tell the program what in fact is water and what is not..
As RSIGNELL explain in viewtopic.php?f=23&t=1438 you can modify your coast data... in my personal experiences I not only need to add some points at the beginning but also at some other "segments" (inserting rows and obligating the coast to follow the trajectory more appropriated for you...) I don't know how to explain it better..sorry...

Good luck: R

[rmueller]

Bug found, captured, and made into butterfly.

I believe the smth_bath call in easygrid should be as follows:
h = smth_bath(h,mask_rho, order, rlim, npass);

instead of
h = smth_bath(h,mask_rho);

The later works, it simply applies the default values every time.

Cheers,
Rachael

[Diego]

Wonderful Rachael!

I'll include this in the next update.

Thanks!

Diego

[rmueller]

Hi Diego,

Another potential contribution is the 'angle' parameter. I see a note in a previous post about adding the scalar value of 'angle' to the output, but no mention of the 'angle' array (unless I missed it elsewhere?). I borrowed some code from the seagrid routine to come up with the following.

Code: Select all

   
% Angles at RHO points
    
lat_u=lat_u.';
lon_u=lon_u.';
[dx,ang]=sw_dist(lat_u(:),lon_u(:),'km');

dx=[dx(:); dx(end)]*1000;  % km==> m
dx=reshape(dx,Lp-1,Mp);
dx=dx.';
dx=[dx(:,1) dx(:,(1:(end-1))) dx(:,end-1)];

ang=[ang(:); ang(end)];
ang=reshape(ang,Lp-1,Mp);
ang=ang.';
ang=[ang(:,1) ang(:,(1:(end-1))) ang(:,end-1)];

    angle = ang*pi/180;   % Radians.

Does this seem right?

Also, aren't x_rho, y_rho, x_psi,... necessary in the netcdf file?

And finally--the answer to this question may be really obvious--what projection are you using? I see that it is not spherical.

Code: Select all

  nc{ 'spherical'}(:)       = 'F';

Do you have suggestions for how to go about applying a different projection (e.g. polar stereographic)? I feel like the projection bit is important, and I still feel quite green in this matter, so your feedback would be helpful!

Thanks,
Rachael

[mistral]

I used this configuration

%--------------------------------------------------------------------------
dims = { 'eta_rho'; 'xi_rho'};
nc{ 'angle'} = ncdouble(dims);
nc{ 'angle'}(:) = angle; % angle of conterclock-wise grid rotation (degrees)
nc{ 'angle'}.description = 'angle of conterclock-wise grid rotation (degrees)';

nc{ 'angle'}.units = 'degrees' ;



this is a good idea???

[jcharris]

Doesn't the section copied below from Utility/get_grid.F indicate that if EASYGRID is modified to store the rotation angle, it should be saved in radians, not degrees?

!
! Read in angle (radians) between XI-axis and EAST at RHO-points.
!
CASE ('angle')
gtype=r2dvar
status=nf_fread2d(ng, model, ncname, ncGRDid(ng), &
& var_name(i), var_id(i), &
& 0, gtype, Vsize, &
& LBi, UBi, LBj, UBj, &
& Fscl, Fmin, Fmax, &
# ifdef MASKING
& GRID(ng) % rmask, &
# endif
& GRID(ng) % angler)

[dnavidad]

hi easygrid-user,

i would like to launch some question about easygrid , particularly with the land process that this program does. Due to i want to use wet_dry option and keep my land i wonder if what it doing easygrid is correct.

i took a look the code to see what is doing with minh, i found that.

Code: Select all

   
h(h<0) = 0;   % Flatten hills and mountains (i.e. positive depths)
h = h + minh; % Add the depth offset minh (specified in USER SETTINGS)
clear xbathy ybathy zbathy
%NOTE: Bathymetry smoothing occurs below "generating masks"

which mean that is taking off any real topography point and doing like a pool in the land with a constant detph equal the value of minh. but that doesn't make sense or that is only an example? it is simple to change it.

1. should not be added this minh to the negative depth preserving the topography? and get the point wet depend on his high and taking into account minh and DCRIT roms parameter.

thanks

[Kirang]

Hello everyone!

Besides EASYGRID and ROMSTOOLS ,is there any other way to create files required to run roms which doesn't require matlab?


Regards,
Katherine

[kate]

Yes, there is a set of python tools which you can download from github. The download command would be:

Code: Select all

git clone https://github.com/kshedstrom/pyroms

[Kirang]

Thank you so much Kate .I will try that.

Regards,
Katherine

[MIDHILAVARNAV]

I tried using make_grid.m of ROMS_tools to create grid, but I did not get any desired result. Kindly resolve the problem.

[johnluick]

GridBuilder - see discussion elsewhere in roms forum.
http://austides.com/downloads/

[andres]

Hi MIDHILAVARNAV,

These tools are used with ROMS_AGRIF/CROCO code, please post your questions here:

https://gforge.inria.fr/forum/forum.php?forum_id=3526

Please post there the crocotools_param.m or romtools_param.m file that you use.