Welcome to mesher’s documentation!

Algorithm

The paper describing mesher can be found here and is the best description for how the algorithm works. However, a brief summary is given here to help the user understand why certain things behave the way they do.

The basic principal of mesher is to generate an unstructured triangular mesh of irregularly shaped and sized triangles. Specifically, as each triangle can be though of as an interpolant to an underlying surface (e.g., DEM), a triangle will therefore approximate the underlying surface to some degree of precision. The main idea of mesher is to generate meshes with fewer computational elements than a raster for use in a hydrological or landsurface model. Importantly, the heterogeneity in key aspects of the surface must be maintained while doing so.

Triangles (starting very large) are inserted into the triangulation, and all the underlying raster cells that are covered by this triangle are used with an error metric, e.g., root mean squared error (RMSE) to determine if this triangle approximate those raster cells “well enough”. “Well enough” is defined by the user as a maximum error tolerance. For example, as user could specify that they wanted a triangle to approximate the underlying raster cells to no more than 5m RMSE. Thus, a triangle is rejected and refined (made smaller) if this tolerance is not fulfilled, i.e., there is too much error introduced. This proceeds for all triangles until all triangles have either a) passed the tolerance or b) the triangle cannot be made smaller, and is accepted.

An example of this iterative refinement of the mesh is shown below for a gaussian hill.

And an example for the Granger sub-basin

As the tolerance is made smaller, less error is introduced, and the mesh better approximates the underlying surface. However, the total number of triangles increases. There is therefore a user’s choice in balancing the total number of triangles, which will require more computational effort to use in a simulation, versus limiting the total error introduced.

Consider the following three meshes where each has a progressively larger tolerance to approximate a base 30 m SRTM DEM.

Therefore, with the tight tolerance of 5 m, a generally uniform distribution of triangle size is created. From a hydrological model’s computational standpoint, having a good mix of triangle sizes is best. This ensure heterogeneity where it is required (smaller triangles), e.g. along ridge lines, and larger triangles where it is less needed, e.g., valleys. In this case, the 5 m tolerance is likely a bit too tight, and will result in limited gains over a raster DEM. The 15 m RMSE tolerance produces a good mix of triangle sizes, where the ridge lines are well represented, but the valleys have few triangles. The 50 m preserves key aspects of the topography, such as the ridge lines, but there is a clear smoothing of the mountains

As mesher is designed for use in hydrological landsurface models, using only the topographic information to create the mesh may be insufficient. Capturing important heterogeneities is key. Therefore, mesher uses a multi-objective approach where it will attempt to fulfill tolerances across multiple input rasters. Shown below is such a case where topography and vegetation height were used. Key aspects of the surface such as the tree line and stream (due to the riparian area vegetation) have smaller triangles to capture the spatial heterogeneity, where as on the gentle and mostly uniform slopes has larger triangles.

Installation

Installation of mesher is possible via pip. Installation into a conda environment probably works but is not tested. Mesher is only supported on Macos and Linux for Python 3.7+

Wheels are not prebuilt for mesher. Instead, mesher will need to be compiled as part of the pip install step. This thus requires a functional build environment.

Mesher is tested on Macos (brew) and Ubuntu (apt-get) although other configurations likely work as expected. Adjust the dependencies below as needed.

It is easiest if Python 3.7, 3.8, 3.9 is used as one of the dependencies (vtk) has prebuilt wheels (see here for details on wheel availability). Consider using the pyenv python version manager if this is not your system default Python version.

Setup environment

Mesher can be built against system libraries or against conan libraries.

Note

Depending on your python install, pip may be pip3

System

Ensure the following are installed via package manager:

For macos:

brew install gdal
brew install boost
brew install cgal
brew install metis

For Ubuntu:

sudo apt-get install libgdal-dev
sudo apt-get install libcgal-dev
sudo apt-get install gdal-bin
sudo apt-get install libcgal-dev
sudo apt-get install libboost-filesystem-dev
sudo apt-get install libboost-program-options-dev
sudo apt-get install libmetis-dev

# on Ubuntu 20.04+
sudo apt-get install python3-gdal
# prior to Ubuntu 20.04, use this instead of python3-gdal
# sudo apt-get install python-gdal

Then install mesher with

pip install mesher

Conan

Install conan via

pip install conan

And then setup a new profile

conan profile new default --detect
conan config install https://github.com/Chrismarsh/conan-config.git

This configuration file setups use of revisions, two new remotes (bincrafters, CHM), and tweaks the settings.yml file to have ubuntu-18.04 and ubuntu-20.04 distros. Setting os.distro = 'ubuntu-20.04' will enable the use of prebuilt library binaries.

Then setup conan to use the new C++ ABI and C++ standard

conan profile update settings.compiler.cppstd=14 default

If using clang (e.g.,Macos), do

conan profile update settings.compiler.libcxx=libc++ default  #with clang

and if using gcc, do

conan profile update settings.compiler.libcxx=libstdc++11 default  #with gcc

then install mesher with

USE_CONAN=TRUE pip install mesher

conda

Warning

This is not tested! Mixing conan + conda seems to not be reliable so please use system libraries.

The Anaconda python environment supports pip installs. This example shows installing Anaconda, however if you already have Anaconda installed, then only the instructions from conda create onward is required.

wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh
bash Miniconda3-latest-Linux-x86_64.sh -b -p $HOME/conda
source $HOME/conda/bin/activate
conda init
conda update -y --all
conda create -y --name mesher python=3.8
conda activate mesher
conda install -c conda-forge p11-kit # reported as required as per https://github.com/Chrismarsh/mesher/issues/20
pip install mesher

This approach will use the system installed gdal.

Install of github branch

You can optionally use pip to install the most recent github version or a github branch. However, the automatic setup of the build environment does not occur, so ensure scikit-build, cmake, conan, and ninja are installed. Then,

pip install git+https://github.com/Chrismarsh/mesher@branch-name

If mesher is already installed, use --force-reinstall to reinstall it.

Compilation

Warning

Using mesher this way (i.e., built from source) will require setting the mesher binary path as described in Environment variables.

Building mesher from source has two parts

compile the backend C++ binary and
setup a python environment for the python frontend.

Setup the build environment as described on the install page

Clone repo

An out of source build is recommended. For example:

cd ~/
git clone https://github.com/Chrismarsh/mesher.
mkdir ~/build && cd ~/build

Python packages

setuptools
wheel
conan
scikit-build>=0.11.1
ninja
vtk
pygdal-chm=="`gdal-config --version`.*"
numpy
scipy
matplotlib
cloudpickle
metis

(Optional) Conan dependencies

Optionally mesher’s dependencies can be build using conan for dependency management.

All of the mesher dependencies are built on Github-CI and uploaded to the bintray repository to serve prebuilt binaries. This means that if the mesher build is done with supported compilers and operating system (described later), the dependencies do not need to be built by the end user.

Warning

The python gdal bindings uses a system-wide gdal rather than the conan gdal the mesher C++ backend links against. This will hopefully be resolved in the future. However, as no data passes between the C++ and Python, having different gdal versions poses no problem.

Warning

Conan and conda don’t seem to consistently work. Use at your own risk.

Setup Conan as described on the installation page.

Install the dependencies into your local conan cache (~/.conan/data)

cd ~/build #if you have not already
conan install ~/mesher -if=. --build missing

The -if=. will produce the FindXXX.cmake files required for the mesher build in the current directory, building missing as needed.

Build

You can set the install prefix to be anywhere, such as shown in the example below

cmake ~/mesher -DCMAKE_INSTALL_PREFIX=/opt/mesher

if you’re using conan,

cmake ~/mesher -DCMAKE_INSTALL_PREFIX=/opt/mesher -DUSE_CONAN=True

then

make install

What to do if things aren’t working

On Macos, homebrew tends to update often. Thus, python packages that are installed against homebrew libraries may break in interesting ways. The most common issue is that the pygdal package hasn’t been update for the newest version of gdal in homebrew. This will manifest as an error like gdal 2.3.4 != gdal 3.2.1 . Unfortunately there is not much that can be until pygdal gets updated. It is best to open an issue either on the mesher or nextgis/pygdal githubs.

If there is an install error about gdal-config missing, please ensure that gdal is installed from your package manager and is available on the path. Running gdal-config --version should produce output

Deployment

Notes for how to deploy to Pypi:

pip install scikit-build
pip install twine
pip install wheel

python setup.py sdist bdist_wheel
twine upload  dist/*

Note that version number needs to be incremented for each Pypi upload

Configuration

Mesher depends heavily upon GDAL to handle the geospatial data and the GDAL python bindings. Mesher’s input rasters can be in any 1-band raster that GDAL can open. Input rasters need to have a coordinate system.

Note

Only the items in the Required section are needed for a valid configuration file. However, most of the power of mesher requires use of the other features.

Required

Configuration files are stored in a python file and passed as an argument to mesher on the command line. For example:

python mesher example_mesher_config.py

Therefore this configuration file must be compliant python code, but as such can contain arbitrary python code. Throughout the configuration file, areas and tolerances are taken to be in metres (m).

A minimum configuration file would be:

dem_filename = '../data/ideal_flat.tif'

max_area= 9999999**2  #Effectively unlimited upper area -- allow tolerance check to refine it further
max_tolerance = 5    #1m max RMSE between triangle and underlying elevation
min_area = 5**2     #triangle area below which we will no longer refine, regardless of max_tolerance
errormetric = 'rmse'
mesher_path = '../mesher'

dem_filename

Type:: path

The base DEM to use. Single band raster. The extent of dem_filename is used to define the area that is meshed. Input parameters are clipped to this extent. However, parameters need not cover the entire extent of the simulation domain. Must have a coordinate system. The coordinate system is used to define the coordinate system of the output mesh. This can be overridden with the wkt_out option – see Coordinate Systems.

This DEM can have missing values. In such a case, the domain to be meshed is determined as the largest continuous area. The use of missing values can be used to produce complex and non-rectangular areas. However, the more complex the domain, the greater the number of triangles required along the edge. See simplify on how to avoid this.

max_area

Type:: double

This is a constraint on the maximum size (m^2) of a triangle. Any triangle that is produce that is above this threshold will be rejected and refined. If you wish to set no upper limit on this, specify some very large number, such max_area = 9999999999999999.

max_tolerance

Type:: double

The maximum difference (vertical distance) between the triangle and the underlying raster as measured by the errormeteric. This can optionally be specified as -1 to skip the tolerance checks, which can be useful in producing uniformly sized triangles.

min_area

Type:: double
Default:: pixel_width * pixel_height

A minimum area (m^2) past which mesher should not refine a triangle further. A good setting is the square area of a DEM cell. This does not mean a triangle won’t be smaller than this; rather, if a triangle is below this threshold it will automatically be accepted as valid. This will override the tolerance setting. For example, if the threshold is 3m^2, and a 2m^2 triangle is checked for validity, it will automatically be accepted, without checking the tolerance. A triangle may end up smaller than this threshold due to other splitting that occurs in order to guarantee triangle quality.

Note

Although there is a default value for this, it really should be specified with a value that makes sense.

errormetric

Type:: string

Assigned an string value that determines the error metric to use.

‘mean_tol’ = Mean elevation difference ‘max_tol’ = Max elevation difference ‘rmse’ = RMSE tolerance

The RMSE produces the best distribution of triangle sizes and does not penalized individual raster cells like max_tol does. RMSE should generally be used.

mesher_path

Type:: string

The mesher script needs to know where the backend mesher executable is located. Optionally use the MESHER_EXE environment variable.

nworkers

Type:: int
Default:: Number of physical CPUs

Mesher will use the number of physical CPUs to parallel various tasks. This allows for fine-tuning that parameter.

nworkers_gdal

Type:: int
Default:: Number of physical CPUs

The initial step where all input parameters and initial conditions are made resampled and clipped to the input DEM extent can be memory heavy. Therefore this limits, separate from the number of workers for the main parameterize loop, how many parallel workers should be used to limit memory consumption.

Environment variables

MESHER_EXE

Type:: string

Instead of specifying a mesher_path in the configuration file, the environment variable MESHER_EXE may be set to the binary. If both MESHER_EXE and mesher_path are defined, the configuration file path takes precedence.

MESHER_NWORKERS

Type:: int

Specifies number of CPUs to use for parallel tasks. This will override a configuration file nworkers value.

MESHER_NWORKERS_GDAL

Type:: int

Specifies number of CPUs to use for parallel gdal tasks. This will override a configuration file nworkers_gdal value.

GDAL_CACHEMAX

Type:: string

GDAL’s raster I/O functionality will, by default, cache reads and writes by up-to 5% of total memory to increase efficiency. When dealing with large rasters and parallel processing, this can result in out-of-memory errors. For example, 32 processes can cache 160% of total available RAM. By default, mesher will internally set the parallel processing cache to be no more than ~30% of total RAM shared across all processes. This is about ~1% per process. For most systems this gives a good blend of performance and memory use. However, if out-of-memory situations occur, try setting this lower or to 0. More information can be found on the gdalwarp documentation and wiki page.

Lloyd iterations

lloyd_itr

Type:: int

Enables n Lloyd iterations. As per the CGAL documentation: the goal of this mesh optimization is to improve the angles inside the mesh, and make them as close as possible to 60 degrees. 100 iterations is a suggested amount. However, please note this does invalidate the numerical guarantees about the minimum amount of error introduced to the mesh.

No Lloyd

Lloyd

verbose

Type:: boolean
Default:: False

Enables verbose output

Outputs

user_output_dir

Type:: path

Sets a user-defined output directory instead of the automatically generated folder name.

reuse_mesh

type:

boolean

default:

False

If a mesh was already generated, and only applying a new parametrization is required, enabling this skips the mesh generation step.

write_shp

Type:: boolean
Default:: True

Writes a .shp file corresponding to the produced triangulation. Is an expensive operation and can slow down production of very large meshes with many parameters.

write_vtu

Type:: boolean
Default:: True

Write a .vtu file that can be view in Paraview.

Coordinate systems

By default mesher uses the coordinate system of the DEM as the output mesh’s coordinate system.

use_input_prj

Type:: boolean
Default:: True

Use the input file’s projection. This is useful for preserving a UTM input.

wkt_out

Type:: string

The output coordinate system can be set using this variable. This needs to be in well known text (wkt) format, e.g.,

wkt_out = "PROJCS[\"North_America_Albers_Equal_Area_Conic\"," \
          "     GEOGCS[\"GCS_North_American_1983\"," \
          "         DATUM[\"North_American_Datum_1983\"," \
          "             SPHEROID[\"GRS_1980\",6378137,298.257222101]]," \
          "         PRIMEM[\"Greenwich\",0]," \
          "         UNIT[\"Degree\",0.017453292519943295]]," \
          "     PROJECTION[\"Albers_Conic_Equal_Area\"]," \
          "     PARAMETER[\"False_Easting\",0]," \
          "     PARAMETER[\"False_Northing\",0]," \
          "     PARAMETER[\"longitude_of_center\",-96]," \
          "     PARAMETER[\"Standard_Parallel_1\",20]," \
          "     PARAMETER[\"Standard_Parallel_2\",60]," \
          "     PARAMETER[\"latitude_of_center\",40]," \
          "     UNIT[\"Meter\",1]," \
          "     AUTHORITY[\"EPSG\",\"102008\"]]"

Domain simplification

simplify

Type:: boolean
Default:: False

As described in the algorithm section, the input DEM defines the area to be meshed. If no-data values are present, then the largest continuous area defines the area to be meshed. As a result, complex basin shapes will likely result in the creation of many triangles along the complex edges. This option can be used to simplify the basin outline. Setting simplify to True requires setting a value for simplify_tol.

Once the domain to be meshed is determined (and is represented by a polyline), this polyline is simplified so-as to have no more than simplify_tol meters of error. By default will enable simplify_buffer. See no_simplify_buffer.

Can only be used with a wkt_out that is a projected CRS.

simplify_tol

Type:: double
Default:: 10 m

The maximum error (m) the polygon simplification of simplify can introduce. Be careful as too high a tolerance will cause triangles to be crated that are outside of the bounds of the raster.

simplify_buffer

Type:: double
Default:: -10 m

Sets a negative buffer (i.e., contracts the meshing domain) to give simplify_tol more room to work. That is, usage of the simplify tolerance without this will likely put triangles outside of the valid data domain. Using this allows for the simplification to result in triangles that exist within the data domain. You can disable the use of this buffer by setting no_simplify_buffer=True. simplify_buffer is enabled by default when simplify=True is given.

no_simplify_buffer

Type:: boolean
Default:: False

Disables simplify_buffer when simplify=True

extent

Type:: list[4]
Default:: []

A large DEM may be subset to [xmin ymin xmax ymax]. These are given in the coordinate system of the input DEM.

clip_to_shp

Type:: str
Default:: None

A large DEM may be clipped to the geometry specified by a shape file. This should be a path to the shape file e.g., `/my/path/geom.shp’. The shp must be georeferenced and will be reprojected by gdalwarp as required.

Domain Repair

fill_holes

Type:: boolean
Default:: False

A 1-pass of the GDAL fill no data algorithm. Maximum search distance is defined as max([pixel_height, pixel_width]) * 5.

Input DEM smoothing

If the DEM quality is poor or if triangles close to the elevation raster cell size is required. If the min_area is approximately equal to the cell size of the raster and tolerance parameter ensures triangles of this size are being produced, then in complex terrain the stair stepping of the raster (due to non-continuous first derivative; i.e., slope) impacts the mesh quality as shown below.

Repeated smoothing iterations can be done. Each smoothing iteration using cubic spline and resamples by iter * smoothing_scaling_factor.

This is the result of 1 smoothing iteration.

This is the result of 2 smoothing iterations.

do_smoothing

Type:: boolean
Default:: False

Smooths the input DEM.

smoothing_scaling_factor

type:

double

default:

2.0

Smoothing factor for above option.

max_smooth_iter

type:

int

default:

1

Number of iterations to smooth over. Each iteration the smoothing magnitude increases by iteration * smoothing_scaling_factor.

Parameters

Parameters, such as vegetation cover, flow accumulation, or soils, can be used in mesher in two ways: 1) have their values assigned to a triangle or 2) also constrain the mesh generation, such as shown with the vegetation heights shown on the algorithm page.

No constraint

Parameters are given by key-value pairs in a dictionary, where value is a dict that contains a file name and an aggregation method:

parameter_files = { 'param_name' : {'file':'file.tif','method':'mean'}}

For example:

parameter_files ={
   'soils' : {'file':'/path/to/soils.tif','method':'mean'},
   'canopy_height' : {'file':'/path/to/veg.tif','method':'mean'}
   }

param_name

Type:: string

This is the name of the parameter as it will appear in the final mesh output. I.e., what you want this to be called for use in a model.

file

Type:: string

This is a fully qualified path to the raster file

method

Type:: string
Default:: ‘mean’

The method controls how the various raster cells that make up the triangle are combined together to give. One of mode or mean (average works too). The mean is the mean of all values, whereas mode takes the value that has the greatest number of cells. The mode is good for classified data, such as soil or vegetation type.

Optionally, method may be a user-specified function that accepts a numpy array and returns one value.

do_cell_resize

Type:: bool
Default:: true

Mesher will resize the cells of a raster to match the input DEM’s resolution. However this can result in producing a lot of small triangles around a coarse raster’s cells if a tolerance for this raster is set. For example, if using a coarse vegetation class with a high resoution DEM, it can produce a lot of edge triangles around the vegetation classes. Therefore, if a tolerance on a parameter layer is set, then do_cell_resize=False is automatically set.

drop

Type:: bool
Default:: false

Some parameters may be used to constrain the mesh during generation but a user may not want it written to file. Setting drop=True will not write this parameter to file.

Classifier

A user-specified function can be given as a classified function. This allows, after the method function has combined all the rasters into a single value for the triangle, to re-value this value. This function takes 1 value, and returns 1 value. It is called after the method function is called.

For example, perhaps a vegetation density metric as derived from remote sensing is to be converted to a canopy height. This could be done as:

def Tree_cover_2_VegHeight(value):
    if value >= 50:
        value = 10 # set the tree hight to 10 m
    else:
        value = 0.1 # otherwise set low density to open
    return value

parameter_files ={'CanopyHeight': {'file':'60N_120W_treecover2010_v3.tif',
                                    'method':'mean',
                                    'classifier':Tree_cover_2_VegHeight
                                    }
                  }

Alternatively, a binary tree/no-tree parameter could be derived

 def Tree_cover_2_Simple_Canopy(value):
     if value >= 50:
         value = 0 # not open
     else:
         value = 1 # open
     return value

parameter_files ={ 'landcover': {'file':'60N_120W_treecover2010_v3.tif',
                                  'method':'mean',
                                  'classifier':Tree_cover_2_Simple_Canopy
                                } }

Multiple input rasters can be combined in to a single parameter using a more complex classifier. This works by passing file and method a list of length n of the files and aggregation methods. The classifier then takes n arguments and returns a singlge value. The following examples shows using a water mask and a landcover map to make a water/open/treed dataset.

def make_landcover(water,veg):
  if water == 1:
    return 0 #keep it water

  if veg < 0.2:
    return 1 #this will be clearing/open

  return 2 # this is tree cover

parameter_files = { 'landcover':{'file':['waterMask.tif','Simard_Pinto_3DGlobalVeg_L3C.tif'],
                                 'method':['mode','mean'],
                                 'classifier':make_landcover}}

Using as constraint

By setting an optional tolerance value in the dictionary, a raster can be used to constrain the mesh generation. If method is mode, then this is a fractional percent of the dominate cell values to cover the triangle area. Otherwise, it is RMSE in the units of the raster’s value.

For example:

parameter_files ={
   'soils' : {'file':'/path/to/soils.tif','method':'mean','tolerance':0.6},
   'canopy_height' : {'file':'/path/to/veg.tif','method':'mean','tolerance':2}
   }

Assuming the soils raster is a classification map, then each triangle must have 60% of one soil type as well as be within 2 m RMSE to the canopy height.

Any number of rasters may have a tolerance. Further, if used with a user-defined classifier, then the tolerance check occurs after the classifier has run

Note

As more tolerances are added, or tolerances become tighter, more and more triangles will be produced. Past a certain point, it does not become meaningful to use an approximating mesh!

Initial conditions

Initial conditions work exactly the same way as paramters. However, they are written to a seperate file than the parameters and are intended to be used as intial conditions in a numerical model. They may be used to constrain the mesh.

initial_conditions = { ... }

Shape file constraints

Shape files may be used to further constrain the mesh, for example to rivers or basin outlines. The line segments in the shape file are treated as barries that triangles cannot cross; thus triangle edges represent the shapefile edges exactly. It may be benificial to simplify these edges somewhat so-as to avoid the creation of many small triangles.

This is further shown in flat_stream.

constraints = { 'river_network' :
               {
                  'file': '../data/Stream.shp'
                  'simplify':1 # will be in shp file's original units
               }
         }

constraints.keyname.simplify

Type:: double

Amount to simplify the shapefile edges by. Measured as maximum error between old and new lines. In the units of the shp file.

MPI

By default Mesher will use MPI to launch the Mesher backend tasks.

MPI_nworkers

Type:: int
Default:: Number of cores on machine (e.g., 4, 8, 10, etc)

Set this to limit the number of processors used.

If Mesher is used on a cluster to process a large domain, the use of a job scheduler, such as SLURM, may be optimal.

MPI_exec_str

Type:: string
Default:: None

Set this to a command to use to invoked the MPI job. For example MPI_exec_str=’./submit_job.sh job.sh’ where submit_job.sh invokes the queue submission, e.g., sbatch “$@” and job.sh contains srun –label –unbuffered python “$@”

The exec string used is f”””{MPI_exec_str} {MPI_do_parameterize_path} pickled_param_args_RANK.pickle False {configfile}””” where MPI_do_parameterize_path holds the path to the helper script that is run with python.

If MPI_exec_str is provided MPI_nworkers must also be provided.

Examples

This describes the validation and example datasets, along with sample mesher inputs.

A zip containing all the example code and data can be found here.

flat

Flat DEM, produces 2 triangles

#Example configuration file using the sample data

dem_filename = '../data/ideal_flat.tif'

max_area= 9999999999999999**2  #Effectively unlimited upper area -- allow tolerance check to refine it further
max_tolerance = 5    # 5m max RMSE between triangle and underlying elevation set to -1 to skip tolerance checks
min_area = 5**2     # triangle area below which we will no longer refine, regardless of max_tolerance

MPI_nworkers = 1

flat_veg

Uses the EOSD dataset to mesh with a mode=0.9 threshold. This captures most of the vegetation patches. ALthough it result in an over-generation of triangles between patches, this is due to producing a good gradation from small to larger triangles.

dem_filename = '../data/ideal_flat.tif'

max_area= 50000**2  #Effectively unlimited upper area -- allow tolerance check to refine it further
max_tolerance = 5    # 5 m max RMSE between triangle and underlying elevation set to -1 to skip tolerance checks
min_area = 30**2     # triangle area below which we will no longer refine, regardless of max_tolerance

parameter_files = {
    'landcover': {'file': '../data/eosd.tif', # vegetation landcover
                  'method': 'mode',
                  'tolerance':.9}
}
MPI_nworkers = 1

gaussian_hill

Using the generated gaussian hill dataset, produces a mesh for a guassian hill.

dem_filename = '../data/ideal_gaussianHill.tif'

max_area= 50000**2  #Effectively unlimited upper area -- allow tolerance check to refine it further
max_tolerance = 5    # 5 m maxe RMSE between triangle and underlying elevation set to -1 to skip tolerance checks
min_area = 30**2     #triangle area below which we will no longer refine, regardless of max_tolerance

MPI_nworkers = 1

ideal_ridge

An idealized ridge line

dem_filename = '../data/ideal_ridge.tif'

max_area= 50000**2  #Effectively unlimited upper area -- allow tolerance check to refine it further
max_tolerance = 5    #5m maxe RMSE between triangle and underlying elevation set to -1 to skip tolerance checks
min_area = 30**2     #triangle area below which we will no longer refine, regardless of max_tolerance

MPI_nworkers = 1

ideal_ridge_low_tol

Same as the ideal_ridge, but with a tighter tolerance. Produces more triangles along the rige to better capture it.

dem_filename = '../data/ideal_ridge.tif'

max_area= 50000**2  #Effectively unlimited upper area -- allow tolerance check to refine it further
max_tolerance = 1    #1m maxe RMSE between triangle and underlying elevation set to -1 to skip tolerance checks
min_area = 30**2     #triangle area below which we will no longer refine, regardless of max_tolerance
MPI_nworkers = 1

uniform

Produces a uniform mesh with area = 100 m x 100 m.

dem_filename = '../data/ideal_flat.tif'

max_area= 9999999999999999**2  #Effectively unlimited upper area -- allow tolerance check to refine it further
max_tolerance = -1    # -1 to skip tolerance checks
min_area = 100**2     #triangle area below which we will no longer refine, regardless of max_tolerance
MPI_nworkers = 1

lloyd

Demonstrates the impact of 100 lloyd optimization iterations on the above uniform domain. Compare to the uniform case

dem_filename = '../data/ideal_flat.tif'

max_area= 9999999999999999**2  #Effectively unlimited upper area -- allow tolerance check to refine it further
max_tolerance = -1    #1 -1 to skip tolerance checks
min_area = 100**2     #triangle area below which we will no longer refine, regardless of max_tolerance

lloyd_itr=100
MPI_nworkers = 1

flat_stream

The flat DEM has been constrained to a stream network input as a shape file. This shows the greater number of triangles near the stream. Simplified stream networks produce fewer triangles along the river constraint.

dem_filename = '../data/ideal_flat.tif'

max_area= 9999999999999999**2  #Effectively unlimited upper area -- allow tolerance check to refine it further
max_tolerance = 5    # 5m max RMSE between triangle and underlying elevation set to -1 to skip tolerance checks
min_area = 5**2     #triangle area below which we will no longer refine, regardless of max_tolerance

constraints = { 'river_network' :
					{
						'file': '../data/Stream.shp',
						'simplify':1 # will be in original projection units
					}
 			}
MPI_nworkers = 1

stream_dem

Same as above, but including the the Granger subset DEM.

dem_filename = '../data/granger1m.tif'

max_area= 99999999**2  #Effectively unlimited upper area -- allow tolerance check to refine it further
max_tolerance = 5    #5 m RMSE
min_area = 25**2     #triangle area below which we will no longer refine, regardless of max_tolerance


constraints = { 'river_network' :
					{
						'file': '../data/Stream.shp'
						# 'simplify':5 # will be in original projection units
					}
 			}

lloyd_itr=1
MPI_nworkers = 1

granger

The Granger subset is used, deriving a mesh from only the elevation map.

dem_filename = '../data/granger1m.tif'

max_area= 99999999**2  #Effectively unlimited upper area -- allow tolerance check to refine it further
max_tolerance = 1    #1 m RMSE
min_area = 5**2     #triangle area below which we will no longer refine, regardless of max_tolerance

lloyd_itr=1
simplify=True
MPI_nworkers = 1

granger_low_veg_weight

The granger subset is used with a low elevation tolerance plus low weights on the vegetation map, showing mesher mostly ignoring the vegetation constraints.

dem_filename = '../data/granger1m.tif'

max_area= 50000**2  #Effectively unlimited upper area -- allow tolerance check to refine it further
max_tolerance = 5    #5m maxe RMSE between triangle and underlying elevation set to -1 to skip tolerance checks
min_area = 30**2     #triangle area below which we will no longer refine, regardless of max_tolerance

use_weights = True
topo_weight=0.8
weight_threshold = 0.8
parameter_files = {
    'landcover': {'file': '../data/eosd.tif',
                  'method': 'mode',
                  'weight':0.2,
                  'tolerance':.9}
}


lloyd_itr=1
simplify=True

granger_high_veg_weight

Same as above but with high weight on the vegetation map, showing the algorithm refining triangles to capture the vegetation patches.

dem_filename = '../data/granger1m.tif'

max_area= 50000**2  #Effectively unlimited upper area -- allow tolerance check to refine it further
max_tolerance = 5    # 5m maxe RMSE between triangle and underlying elevation set to -1 to skip tolerance checks
min_area = 30**2     #triangle area below which we will no longer refine, regardless of max_tolerance

use_weights = True
topo_weight=0.2
weight_threshold = 0.8
parameter_files = {
    'landcover': {'file': '../data/eosd.tif',
                  'method': 'mode',
                  'weight':0.8,
                  'tolerance':.9}
}


lloyd_itr=1
simplify=True
MPI_nworkers = 1

flow_accumulation

This domain shows large-extent meshing with a mountain domain for the Bow Valley region near Canmore, Alberta, Canada. This mesh uses a DEM-derived flow accumulation via RichDEM to ensure the mesh captures the high-accumulation locations. The flow accumulation input is generated with the flow.py script in data folder. This requires the RichDEM Python package to be installed.

dem_filename='../data/chro_extent_lowRes.tif'

max_area=5000**2
max_tolerance=50
min_area=200**2

use_input_prj=False
lloyd_itr=100
simplify=True
simplify_tol=100
simplify_buffer=-50

parameter_files = {    
                     'flow_accumulation':{
                         'file':'../data/flow_accumulation.tif',
                         'method':'mean',
                         'tolerance':50
                     }
                   }


MPI_nworkers = 1

flow_accumulation_granger

Flow accumulation for the smaller Granger subbasin.

dem_filename='../data/granger1m.tif'

max_area=5000**2
max_tolerance=10
min_area=5**2

lloyd_itr=100
simplify=True
simplify_tol=100
simplify_buffer=-50

parameter_files = {    
                     'flow_accumulation':{
                         'file':'../data/flow_accumulation_granger.tif',
                         'method':'mean',
                         'tolerance':500
                     }
                   }


MPI_nworkers = 1

dem_smoothing

A small subset of the above Bow Valley domain is extracted to show the impact of smoothing on the output mesh. If the min_area is approximately equal to the cell size of the raster and tolerance parameter ensures triangles of this size are being produced, then in complex terrain the stair stepping of the raster (due to non-continous first derivative; i.e., slope) impacts the mesh quality as shown below.

mesher has an option to smooth the input DEM to lessen this impact. This is enabled via

do_smoothing = True
max_smooth_iter = 1
smoothing_scaling_factor = 1

Each iteration the smoothing magnitude increases by iteration * smoothing_scaling_factor. This is the result of 1 smoothing iteration.

Subsequent iterations, or increases in smoothing_scaling_factor can reduce the stair-stepping further at the cost of increased smoothing of complex terrain.

do_smoothing = True
max_smooth_iter = 2
smoothing_scaling_factor = 1

dem_filename='../data/chro_small.tif'

max_area=5000**2
max_tolerance=10
min_area=30**2

lloyd_itr = 1
use_input_prj=False
do_smoothing = True
max_smooth_iter = 2
smoothing_scaling_factor = 1

MPI_nworkers = 1

Output

Mesher creates a directory with the same name as the configuration file name or a user-defined name as configured by user_output_dir (see Outputs). Within this folder is a subdirectory with the same name as that of the input DEM.

In the primarily output directory are three JSON format files: - *.ic that contains the initial conditions - *.mesh that contains the primary mesh - *.param that contains the per-triangle parameter values

If no parameter or initial conditions are given, then these files are empty.

In the second directory are the reprojected files (`*_projected`), intermediary files, and the triangulation shape file (`*_USM.shp`). A `*.vtu` file is also created for visualizing in 3D in Paraview.

Note

All indicies are 0-indexed!

.mesh

The primary mesh file’s structure is as follows:

{
  "mesh": {
           "vertex": [ [double,double,double], [double,double,double], ... ],
           "nvertex": int,
           "neigh": [ [int,int,int], [int,int,int], ... ],
           "elem": [ [int,int,int], [int,int,int], ... ],
           },
  "is_geographic": int,
  "proj4": string,
  "UTM_zone": int,
  "nelem": int
}

mesh.vertex

Type:: list of [double,double,double]

List of triangle vertexes. Triangles may share a vertex, thus it isn’t stored twice. A vertex has [x,y,z] coordinates in the output coordinate system.

mesh.nvertex

Type:: int

Number of verticies.

mesh.neigh

Type:: list of [int,int,int]

The neighbours to each triangle. Each int is an index into the elem list.

A value of -1 means “missing”.

mesh.elem

Type:: list of [int,int,int]

Verticies are connected via edges to form triangle elements. These indicies are indexes into the vertex array.

For example, consider the flat mesh. It has 2 triangles in the domain. The .mesh file looks like:

{
    "mesh": {
        "vertex": [
            [
                488479.5,
                6713528.0,
                1000.0
            ],
            [
                490527.5,
                6713528.0,
                1000.0
            ],
            [
                490527.5,
                6711480.0,
                1000.0
            ],
            [
                488479.5,
                6711480.0,
                1000.0
            ]
        ],
        "nvertex": 4,
        "neigh": [
            [
                1,
                -1,
                -1
            ],
            [
                -1,
                0,
                -1
            ]
        ],
        "elem": [
            [
                1,
                0,
                2
            ],
            [
                0,
                3,
                2
            ]
        ],
        "is_geographic": 0,
        "proj4": "+proj=utm +zone=8 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs ",
        "UTM_zone": 8,
        "nelem": 2
    }
}

For example, triangle0 is give by mesh.elem[0]:

This says that to define triangle0, use the 1st, the 0th, and 2nd vertex elements. That is:

triangle0 = [490527.5, 6713528.0, 1000.0], [488479.5, 6713528.0, 1000.0], [490527.5, 6711480.0, 1000.0]

It has a neighbour mesh.neigh[0]

[
   1,
   -1,
   -1
]

which says that edge0 is the only edge with a neighbour, and it is triangle1.

.param

The parameter file’s schema is:

{
    "param1": [
        double/int,
        double/int,,
        ...
    ],
   "param2": [
        double/int,
        double/int,,
        ...
    ]
}

param1

Type:: List of length length(elem)

Each item of the list at index i corresponds to a value for that parameter for triangle_*i*.

For example, that parameter file for ideal is

{
    "area": [
        2097152.0,
        2097152.0
    ]
}

Thus the area of triangle0 is 2097152 m^2 or (1448m)^2. The latter is more easily used for approximate comparisons against raster cell sizes.

.ic

The initial conditions file’s schema is:

{
    "ic1": [
        double/int,
        double/int,,
        ...
    ],
   "ic2": [
        double/int,
        double/int,,
        ...
    ]
}

ic1

Type:: List of length length(elem)

Each item of the list at index i corresponds to a value for that initial condition for triangle_*i*.

.vtu

Within the intermediary folder is a .vtu file. This can be opened in Paraview to provide feedback on how the mesh generation went, and if tweaks should be done.

.shp

Within the intermediary folder is an ESRI .shp file with the _USM.shp suffix. This can be loaded directly into a GIS tool.

Tools

mesh2vtu

Combines and converts a .mesh file and a .param file into a corresponding .vtu file

Usage:

mesh2vtu input.mesh input.param output.vtu

meshmerge

Combines multiple parameter files (.param) into one. All the files are validated to ensure each parameter has the same number of elements.

Usage:

meshmerge [-h] [-o outfile] infile [infile ...]

infile

Type:: string

Path to the n input parameter files to merge into one

outfile

Type:: string

Optional path to the output parameter file. If not given, defaults to a filename of concatenated input files.

meshpermutation

Because of how the triangles are generated and subsequently written to file, triangles that are close in space may not be close in the file. Therefore, if the mesh is used in a numerical model, it may result in inefficient access. This tool computes the permutation of a mesh that minimizes the bandwidth of the connectivity matrix.

This tool adds a new field cell_global_id to the .mesh file. A consuming program can sort the triangles (i.e., mesh.elem see: .mesh) by these ids. Then, spatially close triangles will also be close in memory. This may be visualized by plotting the nearest neighbour connectivity before and after reodering.

Before optimization

After optimization

Usage:

meshpermutation [-h] [-o outfile] [-t type] [-i infile]

outfile

File for output. Overwrite input file

type

rcm Performing RCM bandwidth minimization (default)

nd Performing ND fill-in minimization

metis Performing METIS communication minimization partitioning. Requires -n <value> where <value> is the number of partitions. If this mesh is used later for MPI simulation, this should correspond to the number of MPI ranks that will be used.

n

If metis type is selected, this is required. Number of partitions.

meshstats

Produces a csv file of detailed per-triangle information: the error (‘rmse’ is computed), triangle area (units^2), min and max inner angles (degrees). These values are also written to the _USM.shp file in the mesher output directory.

Expects that the .mesh is in a projected CRS.

Note that the triangles’ RMSE is recomputed, so can be a bit slow on very large meshes.

Note

Requires write_shp=True as the shp file is used to load triangle geometry.

Usage:

meshstats dir

dir

Directory to output directory produced by mesher.

Mesher is a novel multi-objective unstructured mesh generation software that allows mesh generation to be generated from an arbitrary number of hydrologically important features while maintaining a variable spatial resolution. Triangle quality is guaranteed as well as a smooth graduation from small to large triangles. Including these additional features resulted in a better representation of spatial heterogeneity versus classic topography-only mesh generation. The paper describing mesher can be found here

Key points

A novel multi-objective unstructured mesh generation software, Mesher, is presented
Heterogeneity in topography is resolved as well as hydrologically important surface and sub-surface attributes
Spatial heterogeneity is better preserved compared to existing mesh generators

Welcome to mesher’s documentation!

Algorithm

Installation

Setup environment

System

Conan

conda

Install of github branch

Compilation

Clone repo

Python packages

(Optional) Conan dependencies

Build

What to do if things aren’t working

Deployment

Configuration

Required

Environment variables

Lloyd iterations

Outputs

Coordinate systems

Domain simplification

Domain Repair

Input DEM smoothing

Parameters

No constraint

Classifier

Using as constraint

Initial conditions

Shape file constraints

MPI

Examples

flat

flat_veg

gaussian_hill

ideal_ridge

ideal_ridge_low_tol

uniform

lloyd

flat_stream

stream_dem

granger

granger_low_veg_weight

granger_high_veg_weight

flow_accumulation

flow_accumulation_granger

dem_smoothing

Output

.mesh

.param

.ic

.vtu

.shp

Tools

mesh2vtu

meshmerge

meshpermutation

meshstats

Indices and tables