Generating ROMS land mask from GSHHS Global Coastline Datase

[shchepet]

Dear All,

A new automatic tool to construct ROMS land mask from GSHHS
(Global Self-consistent Hierarchical High-resolution Shorelines) dataset
is now available as a compile-once -- use-forever command-line operator.
This is the preferred the way and, in fact superior to the more usual practice
of generating an initial approximation for the mask from bottom topography,
and manually editing it using matlab "editmask" tool -- although some
manual editing may still be needed to resolve "delicate" passages --
99.9% of manual editing can be eliminated. Whoever is interested,
please go ahead and obtain package "tools.tar" from

http://www.atmos.ucla.edu/~alex/ROMS/tools.tar

Make a directory and unpack it there. Refer to "README.compiling" for
help in compiling it. Basically you have to edit makefile header to make sure
that you have correct path for linking to netCDF library (netCDF-4, HDF5-
enabled library is required), you have to create ~/bin directory in your home
directory, and also specify installation path (say /usr/local/bin if you
have write permission in it, but in any case you also need ~/bin as well).
Your may compile the entire package or only tools relevant to land-mask
building: gshhs_to_roms_mask, copymask, single_connect, commit_rmask.


Operator, "gshhs_to_roms_mask", reads GSHHS shoreline data set, generates
"raw" ROMS land mask, and writes it into a separate netCDF file. Usage:

Code: Select all

                gshhs_to_roms_mask roms_grid_file.nc

followed by

Code: Select all

                copymask mask.nc roms_grid_file.nc

followed by some inspection, perhaps editing by "editmask" matlab
tool, and/or

Code: Select all

                single_connect i0 j0 roms_grid_file.nc

where "copymask" and "single_connect" are in this package.

input: (1) "lon_rho","lat_rho" and their dimensions read from
roms_grid_file.nc, and

(2) GSHHS costline data, usually named "gshhs_f.b" from

http://www.ngdc.noaa.gov/mgg/shorelines ... hs/latest/

which can be found inside gshhg-bin-VERSION.zip archive

IMPORTANT: File "gshhs_f.b" must be placed somewhere on
the local computer, after which C-program "read_gshhs.c" (comes
with this package) must be edited to make sure that its line

Code: Select all

   char fname[ ]="/path/to/dir/gshhs_f.b";

points to the correct location. (Here I prefer to hard code the file names,
because there are really no options -- it is just a program to read the specific
file with no alternatives -- coarse resolution coastlines do not count.

output: file "mask.nc" created by this program; it contains an
integer variable "land_rho" with dimensions matching the
above, and land_rho=0 for water points, nonzero-value(s)
(typically +2) for land, but may produce different values
for differend land masses and object, depending on tuning
of this program below.

Here is an example: raw land mask for 1.5km Gulf of Mexico model,
1220 x 898 grid points, straight after gshhs_to_roms_mask roms_grid_file.nc
command

raw_mask.png (9.6 KiB) Viewed 7881 times

The same mask, single_connect -s 500 500 roms_grid_file.nc

single_connetc_mask.png (11.42 KiB) Viewed 7881 times

green areas corresponds to the original land generated by gshhs_to_roms_mask,
yellow-orange areas are the ones blocked by single_connect for the reason that there
is no water passage to these points from the user-specified point i=500, j=500 designated
as "water"; and dark-brown points scattered here and there on the coastline are blocked for
the reason that they are identified as "narrow bays".
Argument -s stands for show, meaning that each point changed from water to land is set
to -1 or -2 instead of 0, so it is possible to visualize what is is doing. All the images presented
here are just screen captures from ncview. Once processed by single_connect, it is guaranteed
that all the blue area is mathematically a single-connected domain. Once you decide that
you are happy with your land mask and do not want to change anything any more, the next
operator commit_rmask grid_file.nc will set all the colored points to "land", and your model
is ready to run.

An enlarged portion of the image above to illustrate which narrow bays are blocked:

single_connetc_mask2.png (7.18 KiB) Viewed 7881 times

...Alexander Shchepetkin

[kate]

Thanks for these! It works great.

May I suggest that Rutgers ROMS compute mask_u/mask_v internally instead of needing to read them?

[shchepet]

Generally you want to keep only what is mathematically necessary in the file,
so not only u_mask, and v_mask should be discarded, but few other items as well.

For a spherical grid this is all what you need:

Code: Select all

netcdf pac16_grid {
dimensions:
        xi_rho = 1330 ;
        xi_u = 1329 ;
        eta_rho = 1218 ;
        eta_v = 1217 ;
variables:
        char spherical ;
                spherical:long_name = "grid type logical switch" ;
                spherical:option_T = "spherical" ;
                spherical:option_F = "cartesian" ;
        double lon_rho(eta_rho, xi_rho) ;
                lon_rho:long_name = "longitude of RHO-points" ;
                lon_rho:units = "degree_east" ;
        double lat_rho(eta_rho, xi_rho) ;
                lat_rho:long_name = "latitude of RHO-points" ;
                lat_rho:units = "degree_north" ;
        double lon_psi(eta_v, xi_u) ;
                lon_psi:long_name = "longitude of PSI-points" ;
                lon_psi:units = "degree_east" ;
        double lat_psi(eta_v, xi_u) ;
                lat_psi:long_name = "latitude of PSI-points" ;
                lat_psi:units = "degree_north" ;
        double pm(eta_rho, xi_rho) ;
                pm:long_name = "curvilinear coordinate metric in XI" ;
                pm:units = "meter-1" ;
        double pn(eta_rho, xi_rho) ;
                pn:long_name = "curvilinear coordinate metric in ETA" ;
                pn:units = "meter-1" ;
        double angle(eta_rho, xi_rho) ;
                angle:long_name = "angle between XI-axis and EAST" ;
                angle:units = "radians" ;
        double f(eta_rho, xi_rho) ;
                f:long_name = "Coriolis parameter at RHO-points" ;
                f:units = "second-1" ;
        double mask_rho(eta_rho, xi_rho) ;
                mask_rho:long_name = "mask on RHO-points" ;
                mask_rho:option_0 = "land" ;
                mask_rho:option_1 = "water" ;
        double hraw(eta_rho, xi_rho) ;
                hraw:long_name = "raw bathymetry at RHO-points" ;
                hraw:units = "meter" ;
        double h(eta_rho, xi_rho) ;
                h:long_name = "bathymetry at RHO-points" ;
                h:units = "meter" ;
}

ever though lon_psi and lat_psi are currently not used by ROMS code or
any supporting software, they may be potentially useful in future.

My code no longer reads Coriolis force f: it simply computes it from lat_rho
(occasionally I compute horizontal component of Coriolis force, hence moving
away from what is known as traditional approximation. So I have to read lat_rho
any way to compute both its sin and cos.

lon_rho is needed only for pedantic way of modeling diurnal cycle in a large
domain: time is universal for the entire grid, but Sun altitude depends on time
and longitude as well. It practice this has very little effect to model results in
comparison just having Sun go up and down at the same time for the entire
domain.

Spherical switch is obviously a legacy: historically PROMS plotting package
based on NGAR Graphics used it to identify spherical grids, but that can be done
by alternatively checking for presence of lon_rho vs. x_rho as coordinate variables.


...as far as gshhs_to_roms_mask.F, check out the effect of its FAST_MODE
CPP-switch inside -- try it with and without it being defined -- it does not
affect the outcome, but affects the speed (note that this switch was introduced
about 10 days ago, but the initial June 24 version of the program did not have it.
It is clear what it does by examining the code inside. This is not a bug fix,
but is rather classical example of cache performance optimization.