Map NetCDF data on OGS surface meshes
This tutorials will show you some examples how to map data from NetCDF file format to OGS surface meshes.
Mapping NetCDF data on
ogs_mesh.vtu can be necessary when you want to use Satellite data or data from other simulation software such as mHM as input data that uses the NetCDF format.
The first example will show you how to map output data from mHM as point data on an OGS mesh.
Example 1: Map mHM output data as point data on OGS mesh
In this example, we will map mHM output data containing soil water content (scw) of Germany onto an
ogs_mesh.vtu of a part of the Selke river catchment using ParaView (
Therefore the NetCDF formatted data has to be imported in ParaView, its coordinate system has to be transformed to the coordinate system of the
ogs_mesh.vtu, then the data has to be mapped (Figure 1).
The corresponding files including a state file for ParaView can be downloaded here, PW: tutorial_swc.
Import NetCDF file:
Let’s get started. At first read in the NetCDF file. Usually NetCDF is a raster based–file format.
This means that a raster (unifrom grid) is defined as an array. Each raster/array cell is defined by coordinates of the cell center and variables (e.g. soil water content) linked to the raster cell.
For time-dependent data, usually one array per time step is defined.
swc.nc. A prompt will asks you to specify a NetCDF reader as there are different NetCDF formats out there. As
swc.nc uses the CF conventions, choose
swc.nc will show up in the pipeline browser. mHM NetCDF files do have unregular cell coordinates, which is not ideal for ParaView to import. We need a little trick.
Spherical coordinates to read in the coordinates as cartesian; click
Apply. Sometimes it may be necessary to first check
Spherical coordinates, click apply, then uncheck
Spherical coordinatesand click
Your imported data should look like in Figure 2, containing longitude and latitude as coordinates.
It’s output type should be
vtkStructuredGrid' . Unfortunately, the cell centered points of the original NetCDF file are now vertices of quadriliteral cells in thevtkStructuredGrid`. The cells are now shifted in x and y direction, so we need another trick later on. Before we need to transfrom the longitude–latitude coordinates.
Transform the coordinate system:
Afterwards, you have to transform the coordiante system from longitude–latitude (EPSG:4326) to Gauss–Kruger (EPSG:5684), which is used by this
ogs_mesh.vtu. Gauss–Kruger–4 is a cartesisn coordinate system that uses m as units. Therfore we will use a
ProgrammableFilter to include a python script. Within the python script we use the library pyproj for the cordinate transformation. You have to install it into your local python folder as there is no option to install in the python version that comes with ParaView—let’s say it did not find any. In case you have a C compiler installed you can install pyproj using pip via
pip install pyproj. Make shure that you have the same python version as used by ParaView. If you do not have any C compiler installed, which can be the case on Windows, then download the corresponding precompiled python wheel file (cpXX in the filename reflects the python version) and install via
pip install pyproj‑2.2.2‑cp38‑cp38‑win_amd64.whl (example). After installing pyproj you gotta tell ParaView were this library is located. Just open a python shell in ParaView and add the directory were you installed pyproj to the system path:
Ok. No right–click on
swc.nc in the pipeline browser and then
Add Filter -> Alphabetical -> ProgrammableFilter.
In the added filter set
Output Data Set Type to
Same as Input add the following script:
from pyproj import Transformer, transform import numpy as np # read in data input = inputs longitude = input.Points[:,0] latitude = input.Points[:,1] # define coordinate transform function by defining the coordinates # systems of input and desired output transf = Transformer.from_crs( "EPSG:4326", "EPSG:5684", always_xy=True) # define x and y as np-arrays x = np.empty(len(longitude), dtype=np.float64) y = np.empty(len(longitude), dtype=np.float64) for i in xrange(0, len(longitude)): dest = transf.transform(longitude[i], latitude[i]) x[i] = dest y[i] = dest # write output output.Points[:,0] = x output.Points[:,1] = y
Important: In python3, the function
xrange was renamed to
range. The output should look as in Figure 3.
Extract grid points: To solve the problem with the shifted cells, we now take the
ProgrammableFilter, extract all points out of it and add new cells based on the extracted points. Just right-click on
ProgrammableFilter in the pipeline browser, then
Add Filter->ExtractSelection. Now select all points on the grid by using the
Select Points on (d) function in the upper panel of the render view and click
Copy Active Selection in the Properties Panel. Then hit
Apply. The output is an
vtkUnstructuredGrid that conists of the same number of points and vertex cells (zero–dimensional cell); there should not be any edges. Okay, now we are ready to import the OGS mesh.
ogs_mesh.vtu: Now you can read in the
File->Open. When viewing the the
ogs_mesh.vtu and the output of
ExtractSelection together, you can see that the OGS mesh is much smaller (Figure 4).
Reduce size of grid to
ogs_mesh.vtu: To save some computation time we will cut out a mesh region that roughly fits the
ogs_mesh.vtu file geometry using the filter
Clip. Just right-click on
ExtractSelection in the pipeline browser, then
Add Filter->Alphabetical->Clip. Set
Origin close to the right side boundary of the
1; 0; 0.
Apply. It wil cut out everything right of the clipping plane. Just repeat this to cut out the remaining sides; hit
Invert when clipping the left and the down side. The output is shown in Figure 4.
Map data onto ogs_mesh.vtu:
a) Finally we can now map the data of the volume on the
ogs_mesh.vtu by adding a
PointDatasetInterpolator filter. Set
Input to the prior filter (either the
Clip or the
Source to the
ogs_mesh.vtu. Use the
VoronoiKernel as Kernel (interpolation algorighm) and set
Null Value to -9999. The
VoronoiKernel uses the value of the closest point in the
Input grid as interpolation value for a point in the
Source grid. Tutorials on other interpolation algorithm will follow.
The output mesh is a
ogs_mesh.vtu with the point data of SWC_L01, SWC_L02 (Figure 5). Note that the mapping also finishes if the both grids do not overlap completely!
b) The new grid does only contain the interpolated point data. To retrieve the original cell and field data of
ogs_mesh.vtu we need to combine both data sets by using the
AppendAttributes filter. Just highlight the original imported
ogs_mesh.vtu and the
PointDatasetInterpolator filter and right click,
Add Filter->AppendAttributes. This should output an
vtkUnstructuredGrid that contains the interpolated point data and all data of the original
To save the file highlight
File->Save Data. By choosing PVD as file format you can export a file series of VTU’s (one VTU per time step) that are linked within a PVD file. By choosing VTU you have the choice to ouput either as file seriess (one VTU per time step) but without a linking PVD or to gather the data of all time steps in one file. The latter option did not work for (ParaView 5.6.2, Ubuntu 18.04 LTS).
Please comment if there are problems or if something is not clear.