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");

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);

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);

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");
}

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();

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");
}

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);
}

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);
}

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);
}

Generated for GDAL by doxygen 1.8.8.