Implementation of GRASS GIS module for 3D raster flow line computation

Title:

Implementation of GRASS GIS module for 3D raster flow line computation

Student:

Anna Petrasova, ​North Carolina State University, Open Source Geospatial Research and Education Laboratory

Organization:

​OSGeo - Open Source Geospatial Foundation

Mentors:

​Helena Mitasova, ​Sören Gebbert

GSoC link:

​abstract

Abstract

GRASS GIS' capabilities of representing and processing data in 3D raster (voxel) format, are rare among other GIS software packages and open opportunities for new applications, especially in research. In this project I would like to extend the 3D raster handling capabilities by developing a module for computing 3D vector flow lines of a voxel representing trivariate continuous field. This new module, r3.flow, will extend the functionality of module r.flow into the third dimension.

Introduction

GRASS GIS' capabilities of representing and processing data in 3D raster (voxel) format, are rare among other GIS software packages and open opportunities for new applications, especially in research. In this project I would like to extend the 3D raster handling capabilities by developing a module for computing 3D vector flow lines of a voxel representing trivariate continuous field. This new module, r3.flow, will extend the functionality of module r.flow into the third dimension.

Background

Although GRASS GIS is mainly used for processing 2-dimensional data, it supports 3 dimensions in both vector and raster data representation. Several 3D raster modules are offered by GRASS GIS and are using its 3D raster library; these modules focus mainly on creating and managing voxel data, computing statistics and converting to different formats. However, there is not many modules which would make use of voxel data and which would help in their analysis. In 3D, we miss standard tools used for topographic analysis such as slope, aspect or flow lines. Today, the importance of such tools increases with the availability of new measurement methods producing multidimensional datasets, such as time-series of soil data, water temperature or salinity. Three-dimensional flow lines can help us study the dynamics of elevation time-series in a form of space-time cube voxel representation. The visualization of such flow lines together with isosurfaces and other geographic data can be done directly in GRASS GIS.

The idea

New module r3.flow will use the flow tracing algorithm based on vector-grid approach implemented in r.flow module. The input aspect maps, represented by horizontal and vertical aspect 3D raster maps will be computed by v.vol.rst module. Alternatively the input aspect maps can be replaced by computation of gradient direction internally. New module r3.flow should offer similar functionality as r.flow does, including upslope/downslope flow lines computation, flow lines density voxel and specifying number of cells between flow lines.

Project plan

date	proposed task
2014-05-19 - 2014-05-23 (week 1)	Write pseudocode for flow lines implementation
2014-05-26 - 2014-05-30 (week 2)	Implementing dummy module which would only read 3D raster and write any vector lines to be able to test results
2014-06-02 - 2014-06-13 (week 3-4)	Implementing the actual algorithm
2014-06-16 - 2014-06-20 (week 05)	More testing and fixing problems
2014-06-23 - 2014-06-27 (week 06)	Implement skip parameter (skip voxels for creating flowlines)
June 23	Mentors and students can begin submitting mid-term evaluations
June 27	Mid-term evaluations deadline
2014-06-30 - 2014-07-04 (week 07)	Implementing upslope/downslope calculation
2014-07-07 - 2014-07-11 (week 08)	Implement flow line density 3D raster
2014-07-14 - 2014-07-18 (week 09)	Implement internal aspect calculation
2014-07-21 - 2014-07-25 (week 10)	More testing and fixing unresolved issues
2014-07-28 - 2014-08-08 (week 11-12)	Depending on situation, either continue with work on r3.flow or start to work on related module r3.slope.aspect which might not be finished within GSoC
2014-08-11 - 2014-08-15 (week 13)	Write code documentation and manual page for r3.flow
August 11	Suggested 'pencils down' date. Take a week to scrub code, write tests, improve documentation, etc.
2014-08-18 - 2014-08-22 (week 14)	Submit evaluation and code to Google
August 18	Firm 'pencils down' date. Mentors, students and organization administrators can begin submitting final evaluations to Google.
August 22	Final evaluation deadline
August 22	Students can begin submitting required code samples to Google

Reports

Week 1

I discussed the implementation with mentors. Based on that I will use similar approach as implemented in VTK, which computes the flowlines by integrating the vector field. To have a first prototype quickly, I will implement it in Python using PyGRASS to write vector lines and Python scripting library (array.py). I will write the Python code in a way which would make it easy to port to C.

I was looking at the VTK implementation and visualized 3D vector field data in Paraview to better understand the algorithm. Basic implementation would look like this:

• input are three 3D rasters (3 components of vector field) and vector with seed points (from where the flowlines are integrated)
• for each point we integrate the flowline based on the velocity field and stopping criterion (time/length)
  • the value of velocity field at a particular point is interpolated from the surrounding cells (voxels) probably using trilinear interpolation (which can be found in vtkVoxel class or on Wikipedia)
  • integration will be done with Runge-Kutta method of 2nd or 4th order.
    • so for example with 2nd order: first compute x(i+1) = x(i) + V(i) * delta_T and then x(i+1) = x(i) + delta_T/2 * (V(i)+ V(i+1)) where V is velocity field, x is location of point on flowline and delta_T is integration step.
• output will be vector of flowlines

I will probably keep the code in sandbox and once it is usable, I will put it in addons

Week 2

I started to implement the prototype in Python:

• trilinear interpolation implemented and tested against SciPy method, unit test is included in the code
• Runge-Kutta integration of 4th order implemented too, but that has to be tested more
• I started r3.flow itself but the prototype is not ready yet.

The code is in sandbox​. Things to consider:

• Soeren suggests using interpolation approach by ​Prickett, 1981 instead of trilinear interpolation from voxel centers.
• r3.flow should be able to compute the gradient so that the input can be also just one 3D raster map instead of 3 maps representing the vector field. For now I could probably use NumPy's method gradient based on central differences. When porting to C, there is ​GPDE library in GRASS.

Week 3

I have a running prototype producing flowlines, I tested it on artificial data created from the first example in r3.gwflow. Here are some results:

The flowlines can be generated with forward or backward integration (upstream or downstream). The integration step can be defined:

• as a time step - it is difficult to guess the right value
• as a lenght step in map units or cells (voxels) - that's now the default. The drawback is that it doesn't behave so well around singular points (where velocity is close to zero) as the time step does.

Another parameter is the maximum number of steps which prevents infinite loops I got around the singular point.

I added another input - besides the vector field components, the input can be only one 3D raster (scalar field) and I compute gradient using NumPy's method gradient.

For the next week, I would like to test it on more complicated data (space-time cube data of terrain evolution) and implement a flow accumulation analogue in 3D. I am still not sure how to proceed, I have to either ensure that the step is so small that I will have at least one point of the flowline in the intersected cell or I have to check between each two steps that I didn't skip any cells on the way.

Also I will look at implementation of Runge Kutta method with adaptive step size which might help with the flowlines near a singular point.