From GRASS-Wiki
Revision as of 18:17, 5 January 2021 by HamishBowman (talk | contribs) (Boundary node attribute preparation)

(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)
Jump to: navigation, search

Interfacing with the ADCIRC coastal ocean model

Grid preparation


OceanMesh2D is an open source (GPL3) generator of 2D triangular finite element meshes suitable for use with ADCIRC that performs an adaptive Delaunay triangulation, iteratively, providing fine mesh detail in areas of numerical complexity. It likes to work in lat/lon coordinate system.

You can output your DEM as NetCDF format with r.out.gdal.

You can create a domain boundary polyline with r.contour. You can then output the domain boundary polyline either as Shapefiles with v.out.ogr, or as two column x,y ASCII data with the v.out.ascii.mat addon module for GRASS 6.

Finding a coastline to use

You can use the r.contour module with your DEM to create a level=0 shoreline.

For those working in areas covering the USA, see the CUSP Coastline wiki page.

Boundary node attribute preparation

You can extract the tidal forcing factor amplitude and phase for boundary nodes needed in the fort.15 setup file from raster maps already loaded into GRASS. In the following example the FES2004 dataset is used. (FES2002 for the K2 constituent as K2 is reputed to be broken in the 2004 edition)

See also the TPXO NetCDF global tide dataset from OSU.

Import FES tidal database data:

# import and split phase into components in=NETCDF:"":spectrum out=fes2004.spectrum -ol

MAPS=`g.mlist rast patt="fes2004*" | sort -n -k 3 -t .`
for MAP in $MAPS ; do
    echo "[$MAP]"
    r.region $MAP s=-90 w=-180
    eval ` -s $MAP`
    g.region rast=$MAP
    g.region n=n-$nsres s=s+$nsres w=w+$ewres
    r.mapcalc "$MAP.crop = $MAP"
    r.region $MAP.crop n=89.9375 s=-89.9375 w=0.0625 e=360.0625
    g.region rast=$MAP.crop
    g.remove $MAP
    g.rename $MAP.crop,$MAP.fixed

spectrum="2N2 K1 K2 M2 M4 Mf Mm Msqm Mtm N2 O1 P1 Q1 S2"
for i in `seq 1 14` ; do
    constit=`echo "$spectrum" | tr ' ' '\n' | head -n $i | tail -n 1`
    g.rename fes2004.Ha.$i.fixed,fes2004.$constit.amplitude
    g.rename fes2004.Hg.$i.fixed,fes2004.$constit.phase

g.rename fes2004.K2.amplitude,fes2004.K2.amplitude.bad
g.rename fes2004.K2.phase,fes2004.K2.phase.bad

for map in `g.mlist rast patt="*.phase" | grep -v bad` ; do
    r.mapcalc << EOF
     $map.u = 1 * cos($map)
     $map.v = 1 * sin($map)

Extract amplitude and phase at open boundary nodes:

constits="M2  N2  S2  K1  O1  P1  K2  2N2  Q1  M4"

for const in $constits ; do
  if [ $const = K2 ] ; then
  echo "[$MAP]"

  v.what.rast vect=adcirc_grid_open_bdy_nodes \
    rast=$MAP.amplitude -pi --q | sort -n > $MAP.ampl.prn

  v.what.rast vect=adcirc_grid_open_bdy_nodes \
    rast=$MAP.phase -p --q | sort -n >  $MAP.phase_raw.prn
  v.what.rast vect=adcirc_grid_open_bdy_nodes \
    rast=$MAP.phase.u -pi --q | sort -n >  $MAP.phase.u.prn
  v.what.rast vect=adcirc_grid_open_bdy_nodes \
    rast=$MAP.phase.v -pi --q | sort -n >  $MAP.phase.v.prn

  echo "% $const: grid open boundary nodes" > $MAP.tidal.prn
  echo "% node ampl phase phase.u phase.v" >> $MAP.tidal.prn
  paste -d'|' $MAP.ampl.prn $MAP.phase_raw.prn \
      $MAP.phase.u.prn $MAP.phase.v.prn | \
    cut -f1,2,4,6,8 -d'|' | \
    tr '|' '\t' >> $MAP.tidal.prn

# check for nulls in the source data.
# Manually replace with values from nearby raster cells using d.what.rast
grep '*' *.tidal.prn

Then, in Matlab or Octave:

M2 = load('fes2004.M2.tidal.prn');
N2 = load('fes2004.N2.tidal.prn');
S2 = load('fes2004.S2.tidal.prn');
K1 = load('fes2004.K1.tidal.prn');
O1 = load('fes2004.O1.tidal.prn');
K2 = load('fes2002.K2.tidal.prn');
P1 = load('fes2004.P1.tidal.prn');
Q1 = load('fes2004.Q1.tidal.prn');

c2N2 = load('fes2004.2N2.tidal.prn');
M4 = load('fes2004.M4.tidal.prn');

% compare nearest neighbor phase angle vs. interpolated phase angle reconstructed from components
[M2(:,3)  cart2pol(M2(:,4), M2(:,5)) * 180/pi]

delta = M2(:,3) - (cart2pol(M2(:,4), M2(:,5)) * 180/pi);
[M2(:,1) delta]
plot(M2(:,1), delta), grid on
xlabel('boundary node')

[M2(:,1) M2(:,3) (cart2pol(M2(:,4), M2(:,5)) * 180/pi) delta]

constits={ 'M2' 'N2' 'S2' 'K1' 'O1' 'P1' 'K2' 'c2N2' 'Q1' 'M4' }

fid = fopen('grid_constits.prn', 'w')

for i = 1: length(constits)
  const = char(constits(i));
  eval(['data = ' const ';'])
  new_phase = cart2pol(data(:,4), data(:,5)) * 180/pi;
  for j = 1:length(new_phase)
    if (new_phase(j) < 0)
      new_phase(j) = new_phase(j) + 360;
  fprintf(fid, ' %s\n', const);
  for j = 1:length(new_phase)
    fprintf(fid, '   %8.6f  %9.5f\n', data(j,2), new_phase(j));


% then manually fix c2N2 -> 2N2

Creating Manning's n nodal attributes

You can load in your fort.14 mesh as points with, then pull it into the Albers location with v.proj, then create a column of n values with 'v.what.rast -i -p' that you can then form into a fort.13 file with UNIX command line power tools such as 'cut' and 'paste'.

Post processing output data

Importing the fort.14 mesh grid

The addon module for GRASS 6 will import a 2.5D triangular mesh from the ADCIRC coastal ocean model into a 3D GRASS vector map.

  • does not yet support loading in node type (boundary, etc.) as point attributes -- database tables can be a bit slow for maps with potentially millions of features. In the mean time here is a little shell script to extract the boundary nodes on the first open boundary into a separate vector points map:
 bdymap=mesh_bdy_pts -p in=$infile out=$inmap
 num_pts=`grep ' = Number of nodes for open boundary' "$infile" | head -n 1 | awk '{print $1}'`
 cat_list=`grep -A"$num_pts" ' = Number of nodes for open boundary' "$infile" | \
    tail -n +2 | tr '\n' ',' | sed -e 's/ //g' -e 's/,$//'`
 v.extract in=$inmap out=$bdymap list="$cat_list"

Importing elevation or velocity reporting stations

With a little awk you can load in the elevation reporting stations for unit 61 with

These shell script commands will take the station positions found in a `elev_stat.151` file and create a text file of points which you can import.

OUT=`basename "$FILE" .151`.dat

echo "# $FILE  `date`" > "$OUT"
head -n1 "$FILE" | sed -e 's/^/#/' >> "$OUT"
tail -n+2 "$FILE" | sed -e 's/^ //' -e 's/ /|/' -e 's/[ ]*! /|/' >> "$OUT" in="$OUT" out=adcirc_stations \
   x=1 y=2 col='lon double, lat double, name varchar(64)'

Visualizing sub-domains

If you wish to have a look at where the sub-domains of a parallel run are, you can use the following Bash script:

mkdir fort14s
for DIR in PE00?? ; do
   cp -a $DIR/fort.14 fort14s/${DIR}.14
cd fort14s
for file in *.14 ; do
   base=`echo "$file" | sed -e 's/PE00//' -e 's/\.14//'`
   echo -n "[$base]" in="$file" out="submesh_$base" --quiet -p in="$file" out="submesh_${base}_pts" --quiet

Here is a loop for displaying 80 sub-meshes, with random dark colors:

for PE in {00..79} ; do
  echo -n "[$PE]"
  R=`expr $RANDOM / 256`
  G=`expr $RANDOM / 256`
  B=`expr $RANDOM / 256`

  d.vect "submesh_$PE" color="$R:$G:$B"

d.vect US_medium_shoreline color=black width=2

d.vect station_map color=black fcolor=red size=8 icon=basic/circle

Importing unit 61, 62, 63, 64 results

For this a useful approach is to import lat/lon coordinates of the nodes (as points) with the AddOn module for GRASS 6, switch over to your target map projection in another GRASS location, pull in the points vector map from the lat/lon location with v.proj, export them with v.out.ascii only keeping the projected points' x and y coordinate columns, tab separated (node order is preserved through all of this), this only needs to be done once.

Next using `ncdump` extract the 'z' or 'u' and 'v' values at each node point, and `paste` those (using the `paste` Unix command line tool) together with the x and y coordinates into a file ready for or Repeat as necessary if the NetCDF output file includes multiple timsteps. All of this can be and should be scripted to save you time.

- addme: bash code + +

d.barb addon module for GRASS 6 for rendering u,v vector arrows