Streamlines in a steady flow field
The analysis of movements in a flow field is a theme of relevance for physics, engineering, Earth and environmental sciences, among others. Two specific examples are ice and atmospheric flows. However no GIS software, to my knowledge, contains modules for pathlines or streamlines calculation.
I present here an algorithm, implement in both Python and C++, that determines pathlines for particles embedded in a 2D, steady, non-uniform flow (example in Fig. 1). This script is independent of any specific GIS software. Its input and output are directly readable by current GIS software (e.g., QuantumGIS, Saga, ArcGis).
Fig 1. Pathlines in the David Glacier (Antarctica) determined by the Python script.
Base map: LANDSAT satellite mosaic. Map created with QGIS Copiaco.
- Language Python 2.6: flowlines_py.zip
- Required modules: Numpy, ogr
- Created and tested with: Python(x,y)
- Version: 2011-10-17
- Language C++: flowlines_cpp.zip
- Required modules (for Windows): gdal14.dll
- Created and tested with: Microsoft Visual C++ 2010 Express
- Version: 2011-10-17
Streamlines, pathlines (and streaklines) are flow representations, that are equivalent in the case of steady flows (Fig. 2). The determination of pathlines can also represent the first step in the determination of timelines, that can be very useful in the analysis of spatial variations of the flows.
Fig. 2. Streamlines (a), pathline (b) e timelines (c) (from Chao-Lung Ting 'Chapter 4: Fluid Kinematics').
What are the methodological requirements to determine the pathlines of a steady vector field?
First, a method to estimate the velocity components for arbitrary points within the field domain. In this script a bilinear interpolation is applied to the components of velocity along the x and y axes.
Second, a method to interpolate the particle movements over time. This is a problem equivalent to solving ordinary differential equations (ODE). One of the most used methods, due to the high-quality results and the possibility of implementing an adaptive step-size, is the 4-5 order Runge-Kutta-Fehlberg method (RKF45).
The 4-5 order Runge-Kutta-Fehlberg method allows to estimate the pathline of particles by using the 5-order RKF algorithm, and the associated errors given a specific time step, by using a 4-order Runge-Kutta estimator. It is thus possible to set a maximum estimation error and to adapt the time step in order to obtain an estimation error lower than the maximum estimation error. This is called adaptive step size control, and is implemented in the current script.
I provide here only the formulas for the order five RKF method (Kreyszig, 2006):
yn+1 = yn + h ( γ1k1 + ... + γ6k6)
with coefficient vector:
γ = [ 16/135 0 6656/12825 28561/56430 -9/50 2/55]
k1 = f(xn, yn)
k2 = f(xn + (1/4) h, yn + h (1/4) k1)
k3 = f(xn + (3/8) h, yn + h ( (3/32) k1 + (9/32) k2) )
k4 = f(xn + (12/13) h, yn + h ( (1932/2197) k1 - (7200/2197) k2 + (7296/2197) k3) )
k5 = f(xn + h, yn + h ( (439/216) k1 - (8) k2 + (3680/513) k3 - (845/4104) k4) )
k6 = f(xn + (1/2) h, yn + h ( - (8/27) k1 + (2) k2 - (3544/2565) k3 + (1859/4104) k4 – (11/40) k5) )
The data input consists of four data types: a parameter files, two grids with vx and vy components and a shapefile with start point for pathlines.
Parameter fileIt is a text file containing the parameters need for the analysis:
- input point shapefile
- input arc/info ascii grid with velocity component along the x axis
- input arc/info ascii grid with velocity component along the y axis
- time step
- total integration time
- tolerance error: numerical value, in the distance unit. Generally a very low number, e.g., 1.0E-6.
- name of output point shapefile
An example of input parameter file is:
- pt_02b.shp // input point shapefile
- vx.asc // ascii grid with x cartesian components of vector field
- vy.asc // ascii grid with y cartesian components of vector field
- 0.01 // time step, in time units
- 100 // total time, in time units
- 1.0e-5 // tolerance error
- pathlines_rkf.shp // output shapefile
Velocity componentsFor the Python implementation, the two velocity component layers must be in Arc/Info ASCII format grids.
In the C++ implementations, all raster format read by GDAL can be used (provided the two grids are in two distinct files).
These two grids must have the same geographical extension and resolution (i.e. coordinates for lower-left cell, cell resolution).
Start pointsOne or more points, stored in a shapefile. Each point must have a unique integer identifier in a field named 'id'.
Point shapefileIt contains the pathlines (as point sequences) for each given start point. Each pathline will be calculated until the total integration time is achieved, or the boundaries of the velocity field are reached. For each calculated point the following information are provided (e.g., Fig. 3):
- id: start point identifier
- x: interpolated point, x coordinate
- y: interpolated point, y coordinate
- vx: interpolated velocity, x component
- vy: interpolated velocity, y component
- vmagn: interpolated velocity, magnitude
- d_time: adaptative time step
- t: time from start
- error: estimated error for interpolation step
Fig. 3. Attribute table for a pathline shapefile.
- Kreyszig, E., 2006. Advanced Engineering Mathematics. John Wiley & Sons.