Programování s knihovnou OGR

Tato stránka obsahuje ukázku zdrojového textu demonstračního programu využívající pro přístup k vektorovým datům knihovnu OGR. Program je inspirován skriptem z cvičení GIS2.

Vypsání atributů

C++

#include <iostream>
#include "ogrsf_frmts.h"

using std::cout;
using std::cerr;
using std::endl;

int main(int argc, char **argv)
{
    const char *filename;
    
    OGRDataSource  *poDS;
    OGRLayer       *poLayer;
    OGRFeatureDefn *poFDefn;

    if (argc != 2) {
	cerr << "Pouziti: " << argv[0] << " shapefile" << endl;
	return 1;
    }
    
    filename = argv[1];
    
    // registrovat dostupne OGR ovladace
    OGRRegisterAll();

    // otevrit ShapeFile pro cteni
    poDS = OGRSFDriverRegistrar::Open(filename, FALSE);
    if (poDS == NULL) {
        cerr << "Otevreni '" << filename << "' selhalo." << endl;
        return 1;
    }
    
    // nacist prvni OGR vrstvu (tj. Shapefile)
    poLayer = poDS->GetLayer(0);
    if (poLayer == NULL) {
        cerr << "Nelze nacist OGR vrstvu." << endl;
        return 1;
    }

    // ziskat informace o vrstve
    poFDefn = poLayer->GetLayerDefn();

    cout << "Detekovana OGR vrstva '" << poFDefn->GetName() << "'.\n\n";

    for(int iField = 0; iField < poFDefn->GetFieldCount(); iField++) {
	// ziskat informace o atributovem sloupci
	OGRFieldDefn *poFieldDefn = poFDefn->GetFieldDefn(iField);
	cout << poFieldDefn->GetNameRef() << endl;
    }

    return 0;
}

Python

import sys
import osgeo.ogr as ogr

def main():
    if len(sys.argv) != 2:
        print >> sys.stderr, "Pouziti: %s shapefile" % sys.argv[0]
        return 1
    
    filename = sys.argv[1]

    # otevrit ShapeFile pro cteni
    ds = ogr.Open(filename)
    if ds is None:
        print >> sys.stderr, "Otevreni '%s' selhalo." % filename
        return 1

    # nacist prvni OGR vrstvu (tj. Shapefile)
    lyr = ds.GetLayer(0)
    if lyr is None:
        print >> sys.stderr, "Nelze nacist OGR vrstvu."
    
    # ziskat informace o vrstve
    feat_defn = lyr.GetLayerDefn()

    print "Detekovana OGR vrstva '%s'.\n" % feat_defn.GetName()

    for i in range(feat_defn.GetFieldCount()):
        # ziskat informace o atributovem sloupci
        field_defn = feat_defn.GetFieldDefn(i)
        print field_defn.GetNameRef()
    
    return 0

if __name__ == "__main__":
    sys.exit(main())

Prostorové predikáty

Pro PostGIS je vhodné nastavit

export PGCLIENTENCODING=WIN1250
export PG_USE_COPY=YES

C++

#include <iostream>
#include <ctime>
#include "ogrsf_frmts.h"
 
using std::cerr;
using std::endl;
using std::cout;
 
OGRDataSource *open_dsn(const char *dsn)
{
    OGRDataSource *poDs;
 
    poDs = OGRSFDriverRegistrar::Open(dsn, FALSE);
    if (!poDs)
        cerr << "Otevreni " << dsn << " selhalo." << std::endl;
 
    return poDs;
}
 
int compute_number(OGRLayer *poLayer1, OGRLayer *poLayer2)
{
    int count, i;
    int nfeat;
    bool found;
 
    OGRFeature *poFeat1, *poFeat2;
    OGRGeometry *poGeom1, *poGeom2;
 
    count = 0;
 
    nfeat = poLayer1->GetFeatureCount();
 
    /* sekvencne cti prvky z vsrtvy 1 (body) */
    poLayer1->ResetReading();
    i = 0;
    while(true) {
        cerr << "\r" << i << "/" << nfeat;
        poFeat1 = poLayer1->GetNextFeature();
        if (!poFeat1)
            break;
        poGeom1 = poFeat1->GetGeometryRef();
        if (!poGeom1 || poGeom1->getGeometryType() != wkbPoint)
            continue;
 
        /* sekvencne cti prvky z vrstvy 2 (polygony) */
        found = false;
        poLayer2->ResetReading();
	poLayer2->SetSpatialFilter(poGeom1);
        while (!found) {
            poFeat2 = poLayer2->GetNextFeature();
            if (!poFeat2)
                break;
            poGeom2 = poFeat2->GetGeometryRef();
            if (!poGeom2 || poGeom2->getGeometryType() != wkbPolygon)
                continue;
 
            if (poGeom1->Within(poGeom2)) {
                found = true;
            }
 
            OGRFeature::DestroyFeature(poFeat2);
        }
        if (found)
            count++;
 
        OGRFeature::DestroyFeature(poFeat1);
        i++;
    }
    cerr << "\r";
 
    return count;
}
 
int main(int argc, char **argv)
{
    const char *dsn1, *dsn2, *relation;
    // const char *pszDriverName = "SQLite";
    const char *pszDriverName = "PostgreSQL";
    char      **papszLCO;
    time_t      start;
    OGRSFDriver    *poDriver;
    OGRDataSource  *poDs1, *poDs2, *poDsSL;
    OGRLayer       *poLayer1, *poLayer2;
 
    if (argc != 3) {
        cerr << "Pouziti: " << argv[0]
             << " <bodova vrstva> <polygonova vrstva>" << endl;
        return 1;
    }
 
    dsn1     = argv[1];
    dsn2     = argv[2];
    relation = argv[3];
 
    // registrovat dostupne OGR ovladace                                                 
    OGRRegisterAll();
 
    cout << "Nacitam data..." << endl;
    // otevrit dsn1, dsn2 pro cteni                                                                           
    poDs1 = open_dsn(dsn1);
    poDs2 = open_dsn(dsn2);
    if (!poDs1 || !poDs2) {
        return 1;
    }
 
    // nacist prvni OGR vrstvu
    poLayer1 = poDs1->GetLayer(0);
    poLayer2 = poDs2->GetLayer(0);
    if (!poLayer1 || !poLayer2 ) {
        cerr << "Nelze nacist OGR vrstvu." << endl;
        return 1;
    }
 
    // konverze dat do databaze SpatiaLite
    cout << "Ukladam data do databaze PostGIS..." << endl;
    poDriver = OGRSFDriverRegistrar::GetRegistrar()->GetDriverByName(pszDriverName);
    if (!poDriver) {
	cerr << "Driver " << pszDriverName << " neni dostupny" << endl;
	return 1;
    }
    // poDsSL = poDriver->Open("tmp.sqlite", TRUE);
    poDsSL = poDriver->Open("PG:dbname=pgis_test", TRUE);
    if (!poDsSL) {
	char **papszDSCO = NULL;
	// papszDSCO = CSLSetNameValue(papszDSCO, "SPATIALITE", "YES");
	// poDsSL = poDriver->CreateDataSource("tmp.sqlite", papszDSCO);
    }
    papszLCO = NULL;
    papszLCO = CSLSetNameValue(papszLCO, "OVERWRITE", "YES");
    poLayer1 = poDsSL->CopyLayer(poLayer1, "points", papszLCO);
    poLayer2 = poDsSL->CopyLayer(poLayer2, "polygons", papszLCO);
 
    cout << "Analyzuji data..." << endl;
    start = time(NULL);
    cout << "Pocet bodu:  " << compute_number(poLayer1, poLayer2)
         << "/" << poLayer1->GetFeatureCount() << endl;
    cout << "Cas vypoctu: " << time(NULL) - start << " sec" << endl;
 
    OGRDataSource::DestroyDataSource(poDsSL);
    OGRDataSource::DestroyDataSource(poDs1);
    OGRDataSource::DestroyDataSource(poDs2);
 
    return 0;
}

Související články

• Programování s knihovnou GDAL
• Ukázka skriptů pro GRASS GIS

Externí odkazy

• OGR API
• OGR API Tutorial
• OGR Architecture
• RFC 8: Developer Guidelines

• Ukázka skriptu pro zpracování katastrálních dat ve formátu VFK