zcobell/adcircmodules

Implementation of routines used to manipulate ADCIRC model files in C++

ADCIRCModules

Implementation of routines used to manipulate ADCIRC model files in C++. The code provides a C++ interface as well as a Python interface.

Python

When the python interface is compiled, you will have an ADCIRCModules python library file, or pyadcircmodules. For Windows, this will by pyadcircmodules.dll, _pyadcircmodules.dylib on Mac systems, and Linux type systems this will be _pyadcircmodules.so. In addition to the library, you will see pyadcircmodules.py. CMake will attempt to install these files to the appropriate directory, however, the user can change this.

import pyadcircmodules
m = pyadcircmodules.Mesh("mymesh.grd")
m.read()

Docker

With the latest versions of the library, there is now a Docker image that can be used for interested parties to sidestep the need to build. One day, there will likely be a python pip package, but for now, this is a good way to remove the entry barrier to use. The Docker image is: zcobell/adcircmodules:v0.5.0-beta.1. It can be retrieved using:

docker pull zcobell/adcircmodules:v0.5.0-beta.1
docker run -it zcobell/adcircmodules:v0.5.0-beta.1

Documentation

There is a Doxygen site available here. The documentation is not yet complete, however, it is slowly progressing. The Python interface is not explicitly documented, however, the Python interaction is generated using SWIG, so the function calls are identical.

Compiling

The build system for this project is CMake

Submodules

Submodules are used in this project. You should initialize the repository using:

git submodule update --init

Note

This project is still under development.

Citation

Please appropriately cite this work in publications, reports, and other source code.

Suggested citation:

@misc{cobell2019adcircmodules,
  title        = {{ADCIRCM}odules: a {C}++ and {P}ython interface for manipulation of {ADCIRC} model data},
  author       = {Cobell, Zachary},
  howpublished = {\url{https://github.com/zcobell/ADCIRCModules}},
  year         = {2020}
}

Credits

Library	Function	Source	Included
Abseil	Abseil is an open-source collection of C++ code (compliant to C++11) designed to augment the C++ standard library.	Visit Website	yes, as submodule
Boost	Used for fast file I/O	Visit Website	yes
cxxopts	Lightweight command line parsing for c++	Visit Website	yes
Date	Howard Hinnant's date library for c++	Visit Website	yes, as submodule
gdal	Used to read raster data to interpolate to meshes	Visit Website	no
HDF5	High-performance data management and storage suite	Visit Website	no
nanoflann	Used for kd-tree searches	Visit Website	yes, as submodule
netCDF	Used to read ADCIRC files formatted with this library	Visit Website	no
proj	PROJ is a generic coordinate transformation software that transforms geospatial coordinates from one coordinate reference system (CRS) to another	Visit Website	yes, but can use external
shapelib	Used to convert meshes to shapefile format	Visit Website	yes
swig	Used to generate an interface to Python	Visit Website	no
indicators	Use for pretty progress bars throughout the code	Visit Website	yes

adcircmodules's Issues

fort13 builder crashing with large landcover datasets

Hey Zach,

I was testing out the fort.13 creator and I was having problems with the code failing on larger datasets (16654 X 12031 pixels). I then tried it with a clipped portion of the larger dataset (3930 X 3784 pixels) and it worked fine. Any thoughts to why this would happen? Does docker have a memory limit? The code starts out fine with the larger data, but then dies at about 15%. Thanks!

cmake newbie problems

Hey Zach,

I have been trying to compile adcircModules for use in python using cmake in Visual Studio 2022 on my desktop and I keep running into errors. I was wondering what your preferred method is to compile this code. Do you use the cmake desktop app directly or go through VS? Thanks for your help!

Johnathan

How to make the directional wind stress reduction nodal attribute

• Firstly, thanks to @zcobell for helping me through this process and writing good, open source code.
• This short write-up is for others who may be wanting to create the directional wind stress reduction nodal attribute that's used inside ADCIRC to take into account roughness elements when forcing your model with certain atmospheric models.
• I realize this is not an "issue" per se but due to the lack of other communication channels I hope this suffices. This could be modified for the wiki site perhaps.

Okay...

• First, you'll need to install ADCIRCmodules. Using CMake this is a breeze. The only difficulty I had was finding OpenMP, which can be used to accelerate this computationally expensive process to calculate the directional wind stress nodal attribute.
• Then you'll need to download the CCAP landcover usage dataset available here. I downloaded the 2016 tif file from there.
• After this, I put together a script that looks something like:

 #include <iostream>
 #include <memory>
 #include "adcircmodules.h"
 
 int main() {
   using namespace Adcirc::Geometry;
   using namespace Adcirc::Interpolation;
   using namespace Adcirc::ModelParameters;
 
   std::unique_ptr<Mesh> m(new Mesh("../binera25m_wFP_tide_only.14"));
   m->read();
   m->defineProjection(4326, true);
 
   std::unique_ptr<NodalAttributes> fort13(new NodalAttributes("../binera25m_wFP_tide_only.13", m.get()));
   fort13->read();
   auto midx2 = fort13->locateAttribute("surface_directional_effective_roughness_length");
 
   std::unique_ptr<Griddata> g(new Griddata(m.get(), "consu_epsg26918.tif", 26918));
 
   for (size_t i = 0; i < m->numNodes(); ++i) {
       //if( i < 50){
         g->setInterpolationFlag(i, Adcirc::Interpolation::Average);
       //}
       //else{
       //  g->setInterpolationFlag(i, Adcirc::Interpolation::NoMethod);
       //}
   }
 
   g->setShowProgressBar(true);
   g->readLookupTable("../CCAP_DWind.table");
   g->setRasterInMemory(false);
   std::vector<std::vector<double>> r = g->computeDirectionalWindReduction(true);
 
   for (size_t i =0; i < m->numNodes(); ++i){
       //if(i < 50){
       //  std::cout<< r[i][0] <<std::endl;
       //}
     fort13->attribute(midx2,i)->setValue(r[i]);
   }
 
   fort13->write("../binera25m_wFP_w_dir_ws.13");
 }

• Note that you'll have to create a dummy fort.13 with the directional wind stress nodal attribute or add on the directional wind stress attribute to an existing fort.13 file.
  Basically, if you're creating a fort.13 from scratch you'll need to have the following in your fort.13 file.

surface_directional_effective_roughness_length 
m
12
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
surface_directional_effective_roughness_length
0

• If you're modifying a fort.13 to support this attribute, then you'll have to add the default values to the top of the file and the user defined values to an appropriate part in the fort.13 file.

• Compile your script (I called it make_dwind.cpp) by linking to the location where the ADCIRCmodule header files are located and link to the shared library.

g++ -I /usr/local/include/adcircmodules/  -ladcircmodules -std=c++14 make_dwind.cpp -o make_dwind.x

• It's important at this point, to stress a couple things.

1. You'll need a look up table (LUT) called something like CCAP_DWind.table which one could save the text below into a file. These values are used to assign roughness values at the nodes so they're important to get right. Place this LUT file in your working directory and make sure the path is correct in your script to-be-compiled.

   2    0.500     :High Intensity Developed
  3    0.390     :Medium Intensity Developed
  4    0.500     :Low Intesnity Developed
  5    0.330     :Developed Open Space
  6    0.060     :Cultivated Land
  7    0.060     :Pasture/Hay
  8    0.040     :Grassland
  9    0.650     :Deciduous Forest
 10    0.720     :Evergreen Forest
 11    0.710     :Mixed Forest
 12    0.120     :Scrub/Shrub
 13    0.550     :Palustrine Forested Wetland
 14    0.110     :Palustrine Scrub/Shrub Wetland
 15    0.110     :Palustrine Emergent Wetland
 16    0.550     :Estuarine Forested Wetland
 17    0.120     :Estuarine Scrub/Schrub Wetland
 18    0.110     :Estuarine Emergent Wetland
 19    0.090     :Unconsolidated Shore
 20    0.090     :Bare Land
 21    0.001     :Open Water
 22    0.040     :Palustrine Aquatic Bed
 23    0.040     :Estuarine Aquatic Bed
*The class has default value(=-9999) will be skiped in mapping

2. Secondly, test the program at first with only a handful of nodes (say 50 or so) by uncommenting the code where I assign the interpolation flag "NoMethod". This way you won't have to wait hours to see a code syntax error appear.
3. Reproject the CCAP tif into a ESPG projection as it has to be in a UTM coordinate system. I chose ESPG:26918 for the NY/NJ area but this depends on where your model is. I accomplished this re-projection of the tif file using gdal's gdalwarp method. Such as:

gdalwarp -t_srs ESPG:21918 source.tif output.tif

After you compile and run the code, you'll see a beautiful progress bar that will give you an estimated time for the process until completion. You can check and visualize this fort13 using OceanMesh2D in MATLAB by first reading the files in.

m = msh('fname','fort.14','aux',{'myfort.13'}); 
plot(m,'type','dir','proj','lamb'); 

Note that by default it plots the infinity norm at each point (the largest directional wind stress reduction at each node of the grid).

which yields something like this:

• A close up
  

check_owi_err: Failed to 'read grid specifications/date in Oceanweather pressure file, domain 1' from './fort.221'

Dear sir,

I prepared pressure and wind files (fort.221 and fort.222) for nws = 30. But I am getting the following error:
Failed to 'read grid specifications/date in Oceanweather pressure file, domain 1' from './fort.221'

lon range = 120:.25:140
lat range = 20:.25:40

grid specifications/date in pressure and wind files:
iLat= 80iLong= 80DX= 0.25DY= 0.25SWLat= 10.25SWLon= 110.25DT= 201209210000

Any suggestions would be appreciated. I added the fort.221 and fort.222 files in the attachment.

fort_221_222.zip

How to make a directional wind reduction attribute?

• Firstly, thanks to @zcobell for helping me through this process and writing good, open source code.
• This short write-up is for others who may be wanting to create the directional wind stress reduction nodal attribute that's used inside ADCIRC to take into account roughness elements when forcing your model with certain atmospheric models.
• I realize this is not an "issue" per se but due to the lack of other communication channels I hope this suffices. This could be modified for the wiki site perhaps.

Okay...

• First, you'll need to install ADCIRCmodules. Using CMake this is a breeze. The only difficulty I had was finding OpenMP, which can be used to accelerate this computationally expensive process to calculate the directional wind stress nodal attribute.
• Then you'll need to download the CCAP landcover usage dataset available here. I downloaded the 2016 tif file from there.
• After this, I put together a script that looks something like:

 #include <iostream>
 #include <memory>
 #include "adcircmodules.h"
 
 int main() {
   using namespace Adcirc::Geometry;
   using namespace Adcirc::Interpolation;
   using namespace Adcirc::ModelParameters;
 
   std::unique_ptr<Mesh> m(new Mesh("../binera25m_wFP_tide_only.14"));
   m->read();
   m->defineProjection(4326, true);
 
   std::unique_ptr<NodalAttributes> fort13(new NodalAttributes("../binera25m_wFP_tide_only.13", m.get()));
   fort13->read();
   auto midx2 = fort13->locateAttribute("surface_directional_effective_roughness_length");
 
   std::unique_ptr<Griddata> g(new Griddata(m.get(), "consu_epsg26918.tif", 26918));
 
   for (size_t i = 0; i < m->numNodes(); ++i) {
       //if( i < 50){
         g->setInterpolationFlag(i, Adcirc::Interpolation::Average);
       //}
       //else{
       //  g->setInterpolationFlag(i, Adcirc::Interpolation::NoMethod);
       //}
   }
 
   g->setShowProgressBar(true);
   g->readLookupTable("../CCAP_DWind.table");
   g->setRasterInMemory(false);
   std::vector<std::vector<double>> r = g->computeDirectionalWindReduction(true);
 
   for (size_t i =0; i < m->numNodes(); ++i){
       //if(i < 50){
       //  std::cout<< r[i][0] <<std::endl;
       //}
     fort13->attribute(midx2,i)->setValue(r[i]);
   }
 
   fort13->write("../binera25m_wFP_w_dir_ws.13");
 }

• Note that you'll have to create a dummy fort.13 with the directional wind stress nodal attribute or add on the directional wind stress attribute to an existing fort.13 file.
  Basically, if you're creating a fort.13 from scratch you'll need to have the following in your fort.13 file.

surface_directional_effective_roughness_length 
m
12
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
surface_directional_effective_roughness_length
0

• If you're modifying a fort.13 to support this attribute, then you'll have to add the default values to the top of the file and the user defined values to an appropriate part in the fort.13 file.

• Compile your script (I called it make_dwind.cpp) by linking to the location where the ADCIRCmodule header files are located and link to the shared library.

g++ -I /usr/local/include/adcircmodules/  -ladcircmodules -std=c++14 make_dwind.cpp -o make_dwind.x

• It's important at this point, to stress a couple things.

1. You'll need a look up table (LUT) called something like CCAP_DWind.table which one could save the text below into a file. These values are used to assign roughness values at the nodes so they're important to get right. Place this LUT file in your working directory and make sure the path is correct in your script to-be-compiled.

   2    0.500     :High Intensity Developed
  3    0.390     :Medium Intensity Developed
  4    0.500     :Low Intesnity Developed
  5    0.330     :Developed Open Space
  6    0.060     :Cultivated Land
  7    0.060     :Pasture/Hay
  8    0.040     :Grassland
  9    0.650     :Deciduous Forest
 10    0.720     :Evergreen Forest
 11    0.710     :Mixed Forest
 12    0.120     :Scrub/Shrub
 13    0.550     :Palustrine Forested Wetland
 14    0.110     :Palustrine Scrub/Shrub Wetland
 15    0.110     :Palustrine Emergent Wetland
 16    0.550     :Estuarine Forested Wetland
 17    0.120     :Estuarine Scrub/Schrub Wetland
 18    0.110     :Estuarine Emergent Wetland
 19    0.090     :Unconsolidated Shore
 20    0.090     :Bare Land
 21    0.001     :Open Water
 22    0.040     :Palustrine Aquatic Bed
 23    0.040     :Estuarine Aquatic Bed
*The class has default value(=-9999) will be skiped in mapping

2. Secondly, test the program at first with only a handful of nodes (say 50 or so) by uncommenting the code where I assign the interpolation flag "NoMethod". This way you won't have to wait hours to see a code syntax error appear.
3. Reproject the CCAP tif into a ESPG projection as it has to be in a UTM coordinate system. I chose ESPG:26918 for the NY/NJ area but this depends on where your model is. I accomplished this re-projection of the tif file using gdal's gdalwarp method. Such as:

gdalwarp -t_srs ESPG:21918 source.tif output.tif

After you compile and run the code, you'll see a beautiful progress bar that will give you an estimated time for the process until completion. You can check and visualize this fort13 using OceanMesh2D in MATLAB by first reading the files in.

m = msh('fname','fort.14','aux',{'myfort.13'}); 
plot(m,'type','dir','proj','lamb'); 

Note that by default it plots the infinity norm at each point (the largest directional wind stress reduction at each node of the grid).

which yields something like this:

• A close up
  

surface canopy coefficient

Hey Zach,

Sorry to keep bring you my problems, but I was wondering what kind of input table I needed to calculate surface canopy coefficient. I tried using the standard computeValuesFromRaster with the DWind conversion table, but it seemed to make everything that had an actual landcover value == 0 in the output. Is there a format for the table I should follow? Like landcover of 10 would be 0 for the canopy coefficient? Is there any example of this conversion table? Thank you for your help!