GDAL
GNM Tutorial

This document is intended to describe using the GNM C++ classes to work with networks. It is advised to read the GNM architecture before to understand the purpose and structure of GNM classes.

Managing networks

In the first example we will create a small water network on the base of the set of spatial data (two shapefiles: pipes and wells which are situated at the GDAL source tree: autotest\gnm\data). The use of the common network format - GNMGdalNetwork class - will allow us to select one of the GDAL-supported vector formats for our network - ESRI Shapefile. After the creation we will build a topology and add some additional data: pumps layer, in order to manually edit network topology.

Initially we register GDAL drivers and create some options (string pairs), which will be passed as parameters during network creation. Here we create a network's name.

#include "gnm.h"
#include <vector>
int main ()
{
char **papszDSCO = NULL;
papszDSCO = CSLAddNameValue(papszDSCO, GNM_MD_NAME, "my_pipes_network");
papszDSCO = CSLAddNameValue(papszDSCO, GNM_MD_SRS, "EPSG:4326");
papszDSCO = CSLAddNameValue(papszDSCO, GNM_MD_DESCR, "My pipes network");
papszDSCO = CSLAddNameValue(papszDSCO, GNM_MD_FORMAT, "ESRI Shapefile");
char ** CSLAddNameValue(char **papszStrList, const char *pszName, const char *pszValue) CPL_WARN_UNUSED_RESULT
Add a new entry to a StringList of "Name=Value" pairs, ("Name:Value" pairs are also supported for bac...
Definition: cpl_string.cpp:1838
void GDALAllRegister(void)
Register all known configured GDAL drivers.
Definition: gdalallregister.cpp:62

Some options are obligatory. The following parameters must be specified during the network creation: the path/name; format of network storage; spatial reference system (EPSG, WKT, etc.). The according dataset with the "network part" will be created and the resulting network will be returned.

GNMGenericNetwork* poDS = (GNMGenericNetwork*) poDriver->Create( "..\\network_data", 0, 0, 0, GDT_Unknown,
papszDSCO );
CSLDestroy(papszDSCO);
GDALDriver * GetDriverByName(const char *)
Fetch a driver based on the short name.
Definition: gdaldrivermanager.cpp:590
Format specific driver.
Definition: gdal_priv.h:1424
GDALDataset * Create(const char *pszName, int nXSize, int nYSize, int nBands, GDALDataType eType, char **papszOptions) CPL_WARN_UNUSED_RESULT
Create a new dataset with this driver.
Definition: gdaldriver.cpp:162
GNM class which represents a geography network of generic format.
Definition: gnm.h:195
void CSLDestroy(char **papszStrList)
Free string list.
Definition: cpl_string.cpp:200
@ GDT_Unknown
Definition: gdal.h:61
GDALDriverManager * GetGDALDriverManager(void)
Fetch the global GDAL driver manager.
Definition: gdaldrivermanager.cpp:97

For now we have a void network consisted of only "system layers". We need to populate it with "class layers" full of features, so we open a certain foreign dataset and copy layers from it to our network. Note, that we use GDALDataset:: methods for working with "class layers", because GNMNetwork inherited from GDALDataset.

GDALDataset *poSrcDS = (GDALDataset*) GDALOpenEx("..\\in_data",
GDAL_OF_VECTOR | GDAL_OF_READONLY, NULL, NULL, NULL );
OGRLayer *poSrcLayer1 = poSrcDS->GetLayerByName("pipes");
OGRLayer *poSrcLayer2 = poSrcDS->GetLayerByName("wells");
poDS->CopyLayer(poSrcLayer1, "pipes");
poDS->CopyLayer(poSrcLayer2, "wells");
GDALClose(poSrcDS);
A set of associated raster bands, usually from one file.
Definition: gdal_priv.h:336
virtual OGRLayer * GetLayerByName(const char *)
Fetch a layer by name.
Definition: gdaldataset.cpp:5179
virtual OGRLayer * CopyLayer(OGRLayer *poSrcLayer, const char *pszNewName, char **papszOptions=nullptr) override
Duplicate an existing layer.
This class represents a layer of simple features, with access methods.
Definition: ogrsf_frmts.h:71
void GDALClose(GDALDatasetH)
Close GDAL dataset.
Definition: gdaldataset.cpp:3588
#define GDAL_OF_READONLY
Open in read-only mode.
Definition: gdal.h:424
#define GDAL_OF_VECTOR
Allow vector drivers to be used.
Definition: gdal.h:448
GDALDatasetH GDALOpenEx(const char *pszFilename, unsigned int nOpenFlags, const char *const *papszAllowedDrivers, const char *const *papszOpenOptions, const char *const *papszSiblingFiles) CPL_WARN_UNUSED_RESULT
Open a raster or vector file as a GDALDataset.
Definition: gdaldataset.cpp:3206

After the successful copying we have the network full of features, but with no topology. The features were added and registered in the network but they are still not connected with each other. Now it is time to build the network topology. There are two ways of doing this in GNM: manually or automatically. In the most cases automatic building is more convenient, while manual is useful for small editings. Automatic building requires some parameters: we must specify which "class layers" will participate in topology building (we select our two layers), a snap tolerance, direct and inverse cost, direction, which is equal 0.00005 in our case. If the building will be successful the network's graph will be filled with the according connections.

printf("\nBuilding network topology ...\n");
char **papszLayers = NULL;
for(int i = 0; i < poDS->GetLayerCount(); ++i)
{
OGRLayer* poLayer = poDS->GetLayer(i);
papszLayers = CSLAddString(papszLayers, poLayer->GetName() );
}
if(poGenericNetwork->ConnectPointsByLines(papszLayers, dfTolerance,
dfDirCost, dfInvCost, eDir) != CE_None )
{
printf("Building topology failed\n");
}
else
{
printf("Topology has been built successfully\n");
}
virtual int GetLayerCount() override
Get the number of layers in this dataset.
virtual OGRLayer * GetLayer(int) override
Fetch a layer by index.
virtual const char * GetName()
Return the layer name.
Definition: ogrlayer.cpp:1728
char ** CSLAddString(char **papszStrList, const char *pszNewString) CPL_WARN_UNUSED_RESULT
Append a string to a StringList and return a pointer to the modified StringList.
Definition: cpl_string.cpp:83

At this point we have a ready network with topological and spatial data, which can be used now for different purposes (analysis, converting into different formats, etc). But sometimes it is necessary to modify some network's data. For example we need to add additional features and attach them to our built topology (modify topology). We create a new "class layer" in the network and add one feature to it.

OGRLayer *poNewLayer = poDS->CreateLayer("pumps", , NULL, wkbPoint, NULL );
if( poNewLayer == NULL )
{
printf( "Layer creation failed.\n" );
exit( 1 );
}
OGRFieldDefn fieldDefn ("pressure",OFTReal);
if( poNewLayer->CreateField( &fieldDefn ) != OGRERR_NONE )
{
printf( "Creating Name field failed.\n" );
exit( 1 );
}
pt.setX(37.291466);
pt.setY(55.828351);
poFeature->SetGeometry(&pt);
if( poNewLayer->CreateFeature( poFeature ) != OGRERR_NONE )
{
printf( "Failed to create feature.\n" );
exit( 1 );
}
GNMGFID gfid = poFeature->GetFID();
virtual OGRLayer * CreateLayer(const char *pszName, OGRSpatialReference *poSpatialRef=nullptr, OGRwkbGeometryType eGType=wkbUnknown, char **papszOptions=nullptr)
This method attempts to create a new layer on the dataset with the indicated name,...
Definition: gdaldataset.cpp:4341
A simple feature, including geometry and attributes.
Definition: ogr_feature.h:355
static OGRFeature * CreateFeature(OGRFeatureDefn *)
Feature factory.
Definition: ogrfeature.cpp:246
static void DestroyFeature(OGRFeature *)
Destroy feature.
Definition: ogrfeature.cpp:282
OGRErr SetGeometry(const OGRGeometry *)
Set feature geometry.
Definition: ogrfeature.cpp:437
GIntBig GetFID() const
Get feature identifier.
Definition: ogr_feature.h:712
Definition of an attribute of an OGRFeatureDefn.
Definition: ogr_feature.h:93
virtual OGRErr CreateField(OGRFieldDefn *poField, int bApproxOK=TRUE)
Create a new field on a layer.
Definition: ogrlayer.cpp:665
OGRErr CreateFeature(OGRFeature *poFeature) CPL_WARN_UNUSED_RESULT
Create and write a new feature within a layer.
Definition: ogrlayer.cpp:626
virtual OGRFeatureDefn * GetLayerDefn()=0
Fetch the schema information for this layer.
Point class.
Definition: ogr_geometry.h:811
void setX(double xIn)
Set x.
Definition: ogr_geometry.h:866
void setY(double yIn)
Set y.
Definition: ogr_geometry.h:870
#define OGRERR_NONE
Success.
Definition: ogr_core.h:292
@ OFTReal
Double Precision floating point.
Definition: ogr_core.h:598
@ wkbPoint
0-dimensional geometric object, standard WKB
Definition: ogr_core.h:321

After the successful creation the feature will be registered in the network and we can connect it with others. There can be two possible ways to do this. In the first case we need a real feature which will be an edge in the connection, while in the second case we do not need such feature, and passing -1 into the ConnectFeatures() method means that the special system edge will be created for this connection and added to the graph automatically. In our case we had added only one point feature and we have not got the line one to be an edge, so we will use the "virtual" connection. We pass the GFID of our point as the source, the GFID of one of the existed features as the target and -1 as the connector. Note that we also set the costs (direct and inverse) and the direction of our edge manually and these values will be written to the graph. When we used the automatic connection (which also uses ConnectFeatures() internally) such vales were set automatically according to the rule which we also set before.

if (poDS->ConnectFeatures(gfid ,63, -1, 5.0, 5.0, GNMDirection_SrcToTgt) != GNMError_None)
{
printf("Can not connect features\n");
}
virtual CPLErr ConnectFeatures(GNMGFID nSrcFID, GNMGFID nTgtFID, GNMGFID nConFID=-1, double dfCost=1, double dfInvCost=1, GNMDirection eDir=GNM_EDGE_DIR_BOTH)
Connects two features via third feature (may be virtual, so the identificator should be -1).

After all we correctly close the network which frees the allocated resources.

GDALClose(poDS);

All in one block:

#include "gnm.h"
#include "gnm_priv.h"
int main ()
{
char **papszDSCO = NULL;
papszDSCO = CSLAddNameValue(papszDSCO, GNM_MD_NAME, "my_pipes_network");
papszDSCO = CSLAddNameValue(papszDSCO, GNM_MD_SRS, "EPSG:4326");
papszDSCO = CSLAddNameValue(papszDSCO, GNM_MD_DESCR, "My pipes network");
papszDSCO = CSLAddNameValue(papszDSCO, GNM_MD_FORMAT, "ESRI Shapefile");
GDALDriver *poDriver = GetGDALDriverManager()->GetDriverByName("GNMFile");
GNMGenericNetwork* poDS = (GNMGenericNetwork*) poDriver->Create( "..\\network_data", 0, 0, 0, GDT_Unknown,
papszDSCO );
CSLDestroy(papszDSCO);
if (poDS == NULL)
{
printf("Failed to create network\n");
exit(1);
}
GDALDataset *poSrcDS = (GDALDataset*) GDALOpenEx("..\\in_data",GDAL_OF_VECTOR | GDAL_OF_READONLY, NULL, NULL, NULL );
if(poSrcDS == NULL)
{
printf("Can not open source dataset at\n");
exit(1);
}
OGRLayer *poSrcLayer1 = poSrcDS->GetLayerByName("pipes");
OGRLayer *poSrcLayer2 = poSrcDS->GetLayerByName("wells");
if (poSrcLayer1 == NULL || poSrcLayer2 == NULL)
{
printf("Can not process layers of source dataset\n");
exit(1);
}
poDS->CopyLayer(poSrcLayer1, "pipes");
poDS->CopyLayer(poSrcLayer2, "wells");
GDALClose(poSrcDS);
printf("\nBuilding network topology ...\n");
char **papszLayers = NULL;
for(int i = 0; i < poDS->GetLayerCount(); ++i)
{
OGRLayer* poLayer = poDS->GetLayer(i);
papszLayers = CSLAddString(papszLayers, poLayer->GetName() );
}
if(poGenericNetwork->ConnectPointsByLines(papszLayers, dfTolerance,
dfDirCost, dfInvCost, eDir) != CE_None )
{
printf("Building topology failed\n");
exit(1);
}
else
{
printf("Topology has been built successfully\n");
}
OGRLayer *poNewLayer = poDS->CreateLayer("pumps", , NULL, wkbPoint, NULL );
if( poNewLayer == NULL )
{
printf( "Layer creation failed.\n" );
exit( 1 );
}
OGRFieldDefn fieldDefn ("pressure",OFTReal);
if( poNewLayer->CreateField( &fieldDefn ) != OGRERR_NONE )
{
printf( "Creating Name field failed.\n" );
exit( 1 );
}
OGRFeature *poFeature = OGRFeature::CreateFeature(poNewLayer->GetLayerDefn());
pt.setX(37.291466);
pt.setY(55.828351);
poFeature->SetGeometry(&pt);
if( poNewLayer->CreateFeature( poFeature ) != OGRERR_NONE )
{
printf( "Failed to create feature.\n" );
exit( 1 );
}
GNMGFID gfid = poFeature->GetFID();
if (poDS->ConnectFeatures(gfid ,63, -1, 5.0, 5.0, GNMDirection_SrcToTgt) != GNMError_None)
{
printf("Can not connect features\n");
}
GDALClose(poDS);
}

Analysing networks

In the second example we will analyse the network which we have built in the first example. We will calculate the shortest path between two points via Dijkstra algorithm performing the feature blockings and saving the resulting path into the file.

Initially we open our network, passing the path to its Shapefile dataset.

#include "gnm.h"
#include "gnm_priv.h"
int main ()
{
GNMGenericNetwork *poNet = (GNMGenericNetwork*) GDALOpenEx("..\\network_data",GDAL_OF_GNM | GDAL_OF_UPDATE, NULL, NULL, NULL );
if(poSrcDS == NULL)
{
printf("Can not open source dataset at\n");
exit(1);
}
#define GDAL_OF_UPDATE
Open in update mode.
Definition: gdal.h:430
#define GDAL_OF_GNM
Allow gnm drivers to be used.
Definition: gdal.h:454

Before any calculations we open the dataset which will hold the layer with the resulting path.

GDALDataset *poResDS;
poResDS = (GDALDataset*) GDALOpenEx("..\\out_data",
NULL, NULL, NULL);
if (poResDS == NULL)
{
printf("Failed to open resulting dataset\n");
exit(1);
}

Finally we use the Dijkstra shortest path method to calculations. This path will be found passing over the blocked feature and saved into internal memory OGRLayer, which we copy to the real dataset. Now it can be visualized by GIS.

OGRLayer *poResLayer = poNet->GetPath(64, 41, GATDijkstraShortestPath, NULL);
if (poResLayer == NULL)
{
printf("Failed to save or calculate path\n");
}
else if (poResDS->CopyLayer(poResLayer, "shp_tutorial.shp") == NULL)
{
printf("Failed to save path to the layer\n");
}
else
{
printf("Path saved successfully\n");
}
GDALClose(poResDS);
poNet->ReleaseResultSet(poRout);
GDALClose(poNet);
}
virtual void ReleaseResultSet(OGRLayer *poResultsSet)
Release results of ExecuteSQL().
Definition: gdaldataset.cpp:6470
virtual OGRLayer * CopyLayer(OGRLayer *poSrcLayer, const char *pszNewName, char **papszOptions=nullptr)
Duplicate an existing layer.
Definition: gdaldataset.cpp:4786
virtual OGRLayer * GetPath(GNMGFID nStartFID, GNMGFID nEndFID, GNMGraphAlgorithmType eAlgorithm, char **papszOptions) override
Create path between start and end GFIDs.

All in one block:

#include "gnm.h"
#include "gnmstdanalysis.h"
int main ()
{
GNMGenericNetwork *poNet = (GNMGenericNetwork*) GDALOpenEx("..\\network_data",
NULL, NULL, NULL );
if(poSrcDS == NULL)
{
printf("Can not open source dataset at\n");
exit(1);
}
GDALDataset *poResDS;
poResDS = (GDALDataset*) GDALOpenEx("..\\out_data",
NULL, NULL, NULL);
if (poResDS == NULL)
{
printf("Failed to open resulting dataset\n");
exit(1);
}
poNet->ChangeBlockState(36, true);
OGRLayer *poResLayer = poNet->GetPath(64, 41, GATDijkstraShortestPath, NULL);
if (poResLayer == NULL)
{
printf("Failed to save or calculate path\n");
}
else if (poResDS->CopyLayer(poResLayer, "shp_tutorial.shp") == NULL)
{
printf("Failed to save path to the layer\n");
}
else
{
printf("Path saved successfully\n");
}
GDALClose(poResDS);
poNet->ReleaseResultSet(poRout);
GDALClose(poNet);
}
virtual CPLErr ChangeBlockState(GNMGFID nFID, bool bIsBlock)
Change the block state of edge or vertex.

Generated for GDAL by doxygen 1.9.4.