/*! \page ogr_apitut OGR API Tutorial This document is intended to document using the OGR C++ classes to read and write data from a file. It is strongly advised that the read first review the OGR Architecture document describing the key classes and their roles in OGR. It also includes code snippets for the corresponding functions in C and Python. \section ogr_apitut_read Reading From OGR For purposes of demonstrating reading with OGR, we will constuct a small utility for dumping point layers from an OGR data source to stdout in comma-delimited format. Initially it is necessary to register all the format drivers that are desired. This is normally accomplished by calling OGRRegisterAll() which registers all format drivers built into GDAL/OGR. In C++ : \code #include "ogrsf_frmts.h" int main() { OGRRegisterAll(); \endcode In C : \code #include "ogr_api.h" int main() { OGRRegisterAll(); \endcode Next we need to open the input OGR datasource. Datasources can be files, RDBMSes, directories full of files, or even remote web services depending on the driver being used. However, the datasource name is always a single string. In this case we are hardcoded to open a particular shapefile. The second argument (FALSE) tells the OGRSFDriverRegistrar::Open() method that we don't require update access. On failure NULL is returned, and we report an error. In C++ : \code OGRDataSource *poDS; poDS = OGRSFDriverRegistrar::Open( "point.shp", FALSE ); if( poDS == NULL ) { printf( "Open failed.\n" ); exit( 1 ); } \endcode In C : \code OGRDataSourceH hDS; hDS = OGROpen( "point.shp", FALSE, NULL ); if( hDS == NULL ) { printf( "Open failed.\n" ); exit( 1 ); } \endcode An OGRDataSource can potentially have many layers associated with it. The number of layers available can be queried with OGRDataSource::GetLayerCount() and individual layers fetched by index using OGRDataSource::GetLayer(). However, we wil just fetch the layer by name. In C++ : \code OGRLayer *poLayer; poLayer = poDS->GetLayerByName( "point" ); \endcode In C : \code OGRLayerH hLayer; hLayer = OGR_DS_GetLayerByName( hDS, "point" ); \endcode Now we want to start reading features from the layer. Before we start we could assign an attribute or spatial filter to the layer to restrict the set of feature we get back, but for now we are interested in getting all features. While it isn't strictly necessary in this circumstance since we are starting fresh with the layer, it is often wise to call OGRLayer::ResetReading() to ensure we are starting at the beginning of the layer. We iterate through all the features in the layer using OGRLayer::GetNextFeature(). It will return NULL when we run out of features. In C++ : \code OGRFeature *poFeature; poLayer->ResetReading(); while( (poFeature = poLayer->GetNextFeature()) != NULL ) { \endcode In C : \code OGRFeatureH hFeature; OGR_L_ResetReading(hLayer); while( (hFeature = OGR_L_GetNextFeature(hLayer)) != NULL ) { \endcode In order to dump all the attribute fields of the feature, it is helpful to get the OGRFeatureDefn. This is an object, associated with the layer, containing the definitions of all the fields. We loop over all the fields, and fetch and report the attributes based on their type. In C++ : \code OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn(); int iField; for( iField = 0; iField < poFDefn->GetFieldCount(); iField++ ) { OGRFieldDefn *poFieldDefn = poFDefn->GetFieldDefn( iField ); if( poFieldDefn->GetType() == OFTInteger ) printf( "%d,", poFeature->GetFieldAsInteger( iField ) ); else if( poFieldDefn->GetType() == OFTReal ) printf( "%.3f,", poFeature->GetFieldAsDouble(iField) ); else if( poFieldDefn->GetType() == OFTString ) printf( "%s,", poFeature->GetFieldAsString(iField) ); else printf( "%s,", poFeature->GetFieldAsString(iField) ); } \endcode In C : \code OGRFeatureDefnH hFDefn = OGR_L_GetLayerDefn(hLayer); int iField; for( iField = 0; iField < OGR_FD_GetFieldCount(hFDefn); iField++ ) { OGRFieldDefnH hFieldDefn = OGR_FD_GetFieldDefn( hFDefn, iField ); if( OGR_Fld_GetType(hFieldDefn) == OFTInteger ) printf( "%d,", OGR_F_GetFieldAsInteger( hFeature, iField ) ); else if( OGR_Fld_GetType(hFieldDefn) == OFTReal ) printf( "%.3f,", OGR_F_GetFieldAsDouble( hFeature, iField) ); else if( OGR_Fld_GetType(hFieldDefn) == OFTString ) printf( "%s,", OGR_F_GetFieldAsString( hFeature, iField) ); else printf( "%s,", OGR_F_GetFieldAsString( hFeature, iField) ); } \endcode There are a few more field types than those explicitly handled above, but a reasonable representation of them can be fetched with the OGRFeature::GetFieldAsString() method. In fact we could shorten the above by using OGRFeature::GetFieldAsString() for all the types. Next we want to extract the geometry from the feature, and write out the point geometry x and y. Geometries are returned as a generic OGRGeometry pointer. We then determine the specific geometry type, and if it is a point, we cast it to point and operate on it. If it is something else we write placeholders. In C++ : \code OGRGeometry *poGeometry; poGeometry = poFeature->GetGeometryRef(); if( poGeometry != NULL && wkbFlatten(poGeometry->getGeometryType()) == wkbPoint ) { OGRPoint *poPoint = (OGRPoint *) poGeometry; printf( "%.3f,%3.f\n", poPoint->getX(), poPoint->getY() ); } else { printf( "no point geometry\n" ); } \endcode In C : \code OGRGeometryH hGeometry; hGeometry = OGR_F_GetGeometryRef(hFeature); if( hGeometry != NULL && wkbFlatten(OGR_G_GetGeometryType(hGeometry)) == wkbPoint ) { printf( "%.3f,%3.f\n", OGR_G_GetX(hGeometry, 0), OGR_G_GetY(hGeometry, 0) ); } else { printf( "no point geometry\n" ); } \endcode The wkbFlatten() macro is used above to convert the type for a wkbPoint25D (a point with a z coordinate) into the base 2D geometry type code (wkbPoint). For each 2D geometry type there is a corresponding 2.5D type code. The 2D and 2.5D geometry cases are handled by the same C++ class, so our code will handle 2D or 3D cases properly. Note that OGRFeature::GetGeometryRef() returns a pointer to the internal geometry owned by the OGRFeature. There we don't actually deleted the return geometry. However, the OGRLayer::GetNextFeature() method returns a copy of the feature that is now owned by us. So at the end of use we must free the feature. We could just "delete" it, but this can cause problems in windows builds where the GDAL DLL has a different "heap" from the main program. To be on the safe side we use a GDAL function to delete the feature. In C++ : \code OGRFeature::DestroyFeature( poFeature ); } \endcode In C : \code OGR_F_Destroy( hFeature ); } \endcode The OGRLayer returned by OGRDataSource::GetLayerByName() is also a reference to an internal layer owned by the OGRDataSource so we don't need to delete it. But we do need to delete the datasource in order to close the input file. Once again we do this with a custom delete method to avoid special win32 heap issus. In C++ : \code OGRDataSource::DestroyDataSource( poDS ); } \endcode In C : \code OGR_DS_Destroy( hDS ); } \endcode All together our program looks like this. In C++ : \code #include "ogrsf_frmts.h" int main() { OGRRegisterAll(); OGRDataSource *poDS; poDS = OGRSFDriverRegistrar::Open( "point.shp", FALSE ); if( poDS == NULL ) { printf( "Open failed.\n" ); exit( 1 ); } OGRLayer *poLayer; poLayer = poDS->GetLayerByName( "point" ); OGRFeature *poFeature; poLayer->ResetReading(); while( (poFeature = poLayer->GetNextFeature()) != NULL ) { OGRFeatureDefn *poFDefn = poLayer->GetLayerDefn(); int iField; for( iField = 0; iField < poFDefn->GetFieldCount(); iField++ ) { OGRFieldDefn *poFieldDefn = poFDefn->GetFieldDefn( iField ); if( poFieldDefn->GetType() == OFTInteger ) printf( "%d,", poFeature->GetFieldAsInteger( iField ) ); else if( poFieldDefn->GetType() == OFTReal ) printf( "%.3f,", poFeature->GetFieldAsDouble(iField) ); else if( poFieldDefn->GetType() == OFTString ) printf( "%s,", poFeature->GetFieldAsString(iField) ); else printf( "%s,", poFeature->GetFieldAsString(iField) ); } OGRGeometry *poGeometry; poGeometry = poFeature->GetGeometryRef(); if( poGeometry != NULL && wkbFlatten(poGeometry->getGeometryType()) == wkbPoint ) { OGRPoint *poPoint = (OGRPoint *) poGeometry; printf( "%.3f,%3.f\n", poPoint->getX(), poPoint->getY() ); } else { printf( "no point geometry\n" ); } OGRFeature::DestroyFeature( poFeature ); } OGRDataSource::DestroyDataSource( poDS ); } \endcode In C : \code #include "ogr_api.h" int main() { OGRRegisterAll(); OGRDataSourceH hDS; OGRLayerH hLayer; OGRFeatureH hFeature; hDS = OGROpen( "point.shp", FALSE, NULL ); if( hDS == NULL ) { printf( "Open failed.\n" ); exit( 1 ); } hLayer = OGR_DS_GetLayerByName( hDS, "point" ); OGR_L_ResetReading(hLayer); while( (hFeature = OGR_L_GetNextFeature(hLayer)) != NULL ) { OGRFeatureDefnH hFDefn; int iField; OGRGeometryH hGeometry; hFDefn = OGR_L_GetLayerDefn(hLayer); for( iField = 0; iField < OGR_FD_GetFieldCount(hFDefn); iField++ ) { OGRFieldDefnH hFieldDefn = OGR_FD_GetFieldDefn( hFDefn, iField ); if( OGR_Fld_GetType(hFieldDefn) == OFTInteger ) printf( "%d,", OGR_F_GetFieldAsInteger( hFeature, iField ) ); else if( OGR_Fld_GetType(hFieldDefn) == OFTReal ) printf( "%.3f,", OGR_F_GetFieldAsDouble( hFeature, iField) ); else if( OGR_Fld_GetType(hFieldDefn) == OFTString ) printf( "%s,", OGR_F_GetFieldAsString( hFeature, iField) ); else printf( "%s,", OGR_F_GetFieldAsString( hFeature, iField) ); } hGeometry = OGR_F_GetGeometryRef(hFeature); if( hGeometry != NULL && wkbFlatten(OGR_G_GetGeometryType(hGeometry)) == wkbPoint ) { printf( "%.3f,%3.f\n", OGR_G_GetX(hGeometry, 0), OGR_G_GetY(hGeometry, 0) ); } else { printf( "no point geometry\n" ); } OGR_F_Destroy( hFeature ); } OGR_DS_Destroy( hDS ); } \endcode In Python: \code import sys import ogr ds = ogr.Open( "point.shp" ) if ds is None: print "Open failed.\n" sys.exit( 1 ) lyr = ds.GetLayerByName( "point" ) lyr.ResetReading() feat = lyr.GetNextFeature() while feat is not None: feat_defn = lyr.GetLayerDefn() for i in range(feat_defn.GetFieldCount()): field_defn = feat_defn.GetFieldDefn(i) # Tests below can be simplified with just : # print feat.GetField(i) if field_defn.GetType() == ogr.OFTInteger: print "%d" % feat.GetFieldAsInteger(i) elif field_defn.GetType() == ogr.OFTReal: print "%.3f" % feat.GetFieldAsDouble(i) elif field_defn.GetType() == ogr.OFTString: print "%s" % feat.GetFieldAsString(i) else: print "%s" % feat.GetFieldAsString(i) geom = feat.GetGeometryRef() if geom is not None and geom.GetGeometryType() == ogr.wkbPoint: print "%.3f, %.3f" % ( geom.GetX(), geom.GetY() ) else: print "no point geometry\n" feat = lyr.GetNextFeature() ds.Destroy() \endcode \section ogr_apitut_write Writing To OGR As an example of writing through OGR, we will do roughly the opposite of the above. A short program that reads comma seperated values from input text will be written to a point shapefile via OGR. As usual, we start by registering all the drivers, and then fetch the Shapefile driver as we will need it to create our output file. In C++ : \code #include "ogrsf_frmts.h" int main() { const char *pszDriverName = "ESRI Shapefile"; OGRSFDriver *poDriver; OGRRegisterAll(); poDriver = OGRSFDriverRegistrar::GetRegistrar()->GetDriverByName( pszDriverName ); if( poDriver == NULL ) { printf( "%s driver not available.\n", pszDriverName ); exit( 1 ); } \endcode In C : \code #include "ogr_api.h" int main() { const char *pszDriverName = "ESRI Shapefile"; OGRSFDriverH hDriver; OGRRegisterAll(); hDriver = OGRGetDriverByName( pszDriverName ); if( hDriver == NULL ) { printf( "%s driver not available.\n", pszDriverName ); exit( 1 ); } \endcode Next we create the datasource. The ESRI Shapefile driver allows us to create a directory full of shapefiles, or a single shapefile as a datasource. In this case we will explicitly create a single file by including the extension in the name. Other drivers behave differently. The second argument to the call is a list of option values, but we will just be using defaults in this case. Details of the options supported are also format specific. In C ++ : \code OGRDataSource *poDS; poDS = poDriver->CreateDataSource( "point_out.shp", NULL ); if( poDS == NULL ) { printf( "Creation of output file failed.\n" ); exit( 1 ); } \endcode In C : \code OGRDataSourceH hDS; hDS = OGR_Dr_CreateDataSource( hDriver, "point_out.shp", NULL ); if( hDS == NULL ) { printf( "Creation of output file failed.\n" ); exit( 1 ); } \endcode Now we create the output layer. In this case since the datasource is a single file, we can only have one layer. We pass wkbPoint to specify the type of geometry supported by this layer. In this case we aren't passing any coordinate system information or other special layer creation options. In C++ : \code OGRLayer *poLayer; poLayer = poDS->CreateLayer( "point_out", NULL, wkbPoint, NULL ); if( poLayer == NULL ) { printf( "Layer creation failed.\n" ); exit( 1 ); } \endcode In C : \code OGRLayerH hLayer; hLayer = OGR_DS_CreateLayer( hDS, "point_out", NULL, wkbPoint, NULL ); if( hLayer == NULL ) { printf( "Layer creation failed.\n" ); exit( 1 ); } \endcode Now that the layer exists, we need to create any attribute fields that should appear on the layer. Fields must be added to the layer before any features are written. To create a field we initialize an OGRField object with the information about the field. In the case of Shapefiles, the field width and precision is significant in the creation of the output .dbf file, so we set it specifically, though generally the defaults are OK. For this example we will just have one attribute, a name string associated with the x,y point. Note that the template OGRField we pass to CreateField() is copied internally. We retain ownership of the object. In C++: \code OGRFieldDefn oField( "Name", OFTString ); oField.SetWidth(32); if( poLayer->CreateField( &oField ) != OGRERR_NONE ) { printf( "Creating Name field failed.\n" ); exit( 1 ); } \endcode In C: \code OGRFieldDefnH hFieldDefn; hFieldDefn = OGR_Fld_Create( "Name", OFTString ); OGR_Fld_SetWidth( hFieldDefn, 32); if( OGR_L_CreateField( hLayer, hFieldDefn, TRUE ) != OGRERR_NONE ) { printf( "Creating Name field failed.\n" ); exit( 1 ); } OGR_Fld_Destroy(hFieldDefn); \endcode The following snipping loops reading lines of the form "x,y,name" from stdin, and parsing them. In C++ and in C : \code double x, y; char szName[33]; while( !feof(stdin) && fscanf( stdin, "%lf,%lf,%32s", &x, &y, szName ) == 3 ) { \endcode To write a feature to disk, we must create a local OGRFeature, set attributes and attach geometry before trying to write it to the layer. It is imperative that this feature be instantiated from the OGRFeatureDefn associated with the layer it will be written to. In C++ : \code OGRFeature *poFeature; poFeature = OGRFeature::CreateFeature( poLayer->GetLayerDefn() ); poFeature->SetField( "Name", szName ); \endcode In C : \code OGRFeatureH hFeature; hFeature = OGR_F_Create( OGR_L_GetLayerDefn( hLayer ) ); OGR_F_SetFieldString( hFeature, OGR_F_GetFieldIndex(hFeature, "Name"), szName ); \endcode We create a local geometry object, and assign its copy (indirectly) to the feature. The OGRFeature::SetGeometryDirectly() differs from OGRFeature::SetGeometry() in that the direct method gives ownership of the geometry to the feature. This is generally more efficient as it avoids an extra deep object copy of the geometry. In C++ : \code OGRPoint pt; pt.setX( x ); pt.setY( y ); poFeature->SetGeometry( &pt ); \endcode In C : \code OGRGeometryH hPt; hPt = OGR_G_CreateGeometry(wkbPoint); OGR_G_SetPoint_2D(hPt, 0, x, y); OGR_F_SetGeometry( hFeature, hPt ); OGR_G_DestroyGeometry(hPt); \endcode Now we create a feature in the file. The OGRLayer::CreateFeature() does not take ownership of our feature so we clean it up when done with it. In C++ : \code if( poLayer->CreateFeature( poFeature ) != OGRERR_NONE ) { printf( "Failed to create feature in shapefile.\n" ); exit( 1 ); } OGRFeature::DestroyFeature( poFeature ); } \endcode In C : \code if( OGR_L_CreateFeature( hLayer, hFeature ) != OGRERR_NONE ) { printf( "Failed to create feature in shapefile.\n" ); exit( 1 ); } OGR_F_Destroy( hFeature ); } \endcode Finally we need to close down the datasource in order to ensure headers are written out in an orderly way and all resources are recovered. In C++ : \code OGRDataSource::DestroyDataSource( poDS ); } \endcode In C : \code OGR_DS_Destroy( hDS ); } \endcode The same program all in one block looks like this: In C++ : \code #include "ogrsf_frmts.h" int main() { const char *pszDriverName = "ESRI Shapefile"; OGRSFDriver *poDriver; OGRRegisterAll(); poDriver = OGRSFDriverRegistrar::GetRegistrar()->GetDriverByName( pszDriverName ); if( poDriver == NULL ) { printf( "%s driver not available.\n", pszDriverName ); exit( 1 ); } OGRDataSource *poDS; poDS = poDriver->CreateDataSource( "point_out.shp", NULL ); if( poDS == NULL ) { printf( "Creation of output file failed.\n" ); exit( 1 ); } OGRLayer *poLayer; poLayer = poDS->CreateLayer( "point_out", NULL, wkbPoint, NULL ); if( poLayer == NULL ) { printf( "Layer creation failed.\n" ); exit( 1 ); } OGRFieldDefn oField( "Name", OFTString ); oField.SetWidth(32); if( poLayer->CreateField( &oField ) != OGRERR_NONE ) { printf( "Creating Name field failed.\n" ); exit( 1 ); } double x, y; char szName[33]; while( !feof(stdin) && fscanf( stdin, "%lf,%lf,%32s", &x, &y, szName ) == 3 ) { OGRFeature *poFeature; poFeature = OGRFeature::CreateFeature( poLayer->GetLayerDefn() ); poFeature->SetField( "Name", szName ); OGRPoint pt; pt.setX( x ); pt.setY( y ); poFeature->SetGeometry( &pt ); if( poLayer->CreateFeature( poFeature ) != OGRERR_NONE ) { printf( "Failed to create feature in shapefile.\n" ); exit( 1 ); } OGRFeature::DestroyFeature( poFeature ); } OGRDataSource::DestroyDataSource( poDS ); } \endcode In C : \code #include "ogr_api.h" int main() { const char *pszDriverName = "ESRI Shapefile"; OGRSFDriverH hDriver; OGRDataSourceH hDS; OGRLayerH hLayer; OGRFieldDefnH hFieldDefn; double x, y; char szName[33]; OGRRegisterAll(); hDriver = OGRGetDriverByName( pszDriverName ); if( hDriver == NULL ) { printf( "%s driver not available.\n", pszDriverName ); exit( 1 ); } hDS = OGR_Dr_CreateDataSource( hDriver, "point_out.shp", NULL ); if( hDS == NULL ) { printf( "Creation of output file failed.\n" ); exit( 1 ); } hLayer = OGR_DS_CreateLayer( hDS, "point_out", NULL, wkbPoint, NULL ); if( hLayer == NULL ) { printf( "Layer creation failed.\n" ); exit( 1 ); } hFieldDefn = OGR_Fld_Create( "Name", OFTString ); OGR_Fld_SetWidth( hFieldDefn, 32); if( OGR_L_CreateField( hLayer, hFieldDefn, TRUE ) != OGRERR_NONE ) { printf( "Creating Name field failed.\n" ); exit( 1 ); } OGR_Fld_Destroy(hFieldDefn); while( !feof(stdin) && fscanf( stdin, "%lf,%lf,%32s", &x, &y, szName ) == 3 ) { OGRFeatureH hFeature; OGRGeometryH hPt; hFeature = OGR_F_Create( OGR_L_GetLayerDefn( hLayer ) ); OGR_F_SetFieldString( hFeature, OGR_F_GetFieldIndex(hFeature, "Name"), szName ); hPt = OGR_G_CreateGeometry(wkbPoint); OGR_G_SetPoint_2D(hPt, 0, x, y); OGR_F_SetGeometry( hFeature, hPt ); OGR_G_DestroyGeometry(hPt); if( OGR_L_CreateFeature( hLayer, hFeature ) != OGRERR_NONE ) { printf( "Failed to create feature in shapefile.\n" ); exit( 1 ); } OGR_F_Destroy( hFeature ); } OGR_DS_Destroy( hDS ); } \endcode In Python : \code import sys import ogr import string driverName = "ESRI Shapefile" drv = ogr.GetDriverByName( driverName ) if drv is None: print "%s driver not available.\n" % driverName sys.exit( 1 ) ds = drv.CreateDataSource( "point_out.shp" ) if ds is None: print "Creation of output file failed.\n" sys.exit( 1 ) lyr = ds.CreateLayer( "point_out", None, ogr.wkbPoint ) if lyr is None: print "Layer creation failed.\n" sys.exit( 1 ) field_defn = ogr.FieldDefn( "Name", ogr.OFTString ) field_defn.SetWidth( 32 ) if lyr.CreateField ( field_defn ) != 0: print "Creating Name field failed.\n" sys.exit( 1 ) # Expected format of user input: x y name linestring = raw_input() linelist = string.split(linestring) while len(linelist) == 3: x = float(linelist[0]) y = float(linelist[1]) name = linelist[2] feat = ogr.Feature( lyr.GetLayerDefn()) feat.SetField( "Name", name ) pt = ogr.Geometry(ogr.wkbPoint) pt.SetPoint_2D(0, x, y) feat.SetGeometry(pt) if lyr.CreateFeature(feat) != 0: print "Failed to create feature in shapefile.\n" sys.exit( 1 ) feat.Destroy() linestring = raw_input() linelist = string.split(linestring) ds.Destroy() \endcode */