Programování s knihovnou GDAL

Tato stránka obsahuje ukázku zdrojového textu demonstračního programu využívající pro přístup k rastrovým datům knihovnu GDAL. Program počítá normalizovaný vegetační index, více informace na cvičení předmětu Zpracování obrazových dat.

$ndvi={\frac {(tm4-tm3)}{(tm4+tm3)}}$

kde tm4 je 4. kanál družicové scény Landsat TM, tm3 je kanál třetí.

C++

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

using std::cout;
using std::cerr;
using std::endl;
 
GDALDataset *open_file(const char *filename)
{
    GDALDataset *poDs;
    
    poDs = (GDALDataset *) GDALOpen(filename, GA_ReadOnly);
    if (poDs == NULL) {
        cerr << "Otevreni '" << filename << "' selhalo." << endl;
        return NULL;
    }
    
    return poDs;
}

int calculate_ndvi(GDALDataset *poDsNdvi, GDALDataset *poDsTm3, GDALDataset *poDsTm4)
{
    GDALRasterBand *poBandNdvi, *poBandTm3, *poBandTm4;
    
    poBandNdvi = poDsNdvi->GetRasterBand(1);
    poBandTm3  = poDsTm3->GetRasterBand(1);
    poBandTm4  = poDsTm4->GetRasterBand(1);

    if (!poBandNdvi || !poBandTm3 || !poBandTm4)
	return -1;

    unsigned char *pafScanlineTm3, *pafScanlineTm4;
    double *pafScanlineNdvi;
    int nXSize = poBandTm3->GetXSize();
    int nYSize = poBandTm3->GetYSize();
    if (nXSize != poBandTm4->GetXSize() ||
	nYSize != poBandTm4->GetYSize())
	return -1;

    pafScanlineTm3 = (unsigned char *) CPLMalloc(sizeof(unsigned char) * nXSize);
    pafScanlineTm4 = (unsigned char *) CPLMalloc(sizeof(unsigned char) * nXSize);
    pafScanlineNdvi = (double *) CPLMalloc(sizeof(double) * nXSize);

    for (int row = 0; row < nYSize; row++) {
	poBandTm3->RasterIO(GF_Read, 0, row, nXSize, 1, 
			    pafScanlineTm3, nXSize, 1, GDT_Byte, 
			    0, 0);
	poBandTm4->RasterIO(GF_Read, 0, row, nXSize, 1, 
			    pafScanlineTm4, nXSize, 1, GDT_Byte, 
			    0, 0);
	
	for (int col = 0; col < nXSize; col++) {
	    pafScanlineNdvi[col] = double(pafScanlineTm4[col] - pafScanlineTm3[col]) /
		(pafScanlineTm4[col] + pafScanlineTm3[col]);
	}
	poBandNdvi->RasterIO(GF_Write, 0, row, nXSize, 1, 
			     pafScanlineNdvi, nXSize, 1, GDT_Float32, 
			     0, 0);
    }
    
    CPLFree(pafScanlineTm3);
    CPLFree(pafScanlineTm4);
    CPLFree(pafScanlineNdvi);
    return 0;
}

int main(int argc, char **argv)
{
    const char *filenameTm3, *filenameTm4;
     
    double trans[6];
    GDALDriver   *poDriver;
    GDALDataset  *poDsTm3, *poDsTm4, *poDsNdvi;
    
    if (argc != 3) {
	cerr << "Pouziti: " << argv[0] << " tm3 tm4" << endl;
	return 1;
    }
    
    filenameTm3 = argv[1];
    filenameTm4 = argv[2];

    // registrovat dostupne GDAL ovladace
    GDALAllRegister();
     
    // otevrit rastrove soubory 'tm3' a 'tm4' pro cteni
    poDsTm3 = open_file(filenameTm3);
    poDsTm4 = open_file(filenameTm4);
    if (poDsTm3 == NULL || poDsTm4 == NULL)
	return 1;

    poDriver = poDsTm3->GetDriver();
    poDsNdvi = poDriver->Create("ndvi.tif", poDsTm3->GetRasterXSize(), poDsTm3->GetRasterYSize(),
				1, GDT_Float32, NULL);
    poDsTm3->GetGeoTransform(trans);
    poDsNdvi->SetGeoTransform(trans);
    poDsNdvi->SetProjection(poDsTm3->GetProjectionRef());
    
    if (poDsNdvi == NULL)
	cerr << "Nelze vytvorit vystupni soubor 'ndvi.tif'" << endl;

    if (calculate_ndvi(poDsNdvi, poDsTm3, poDsTm4) != 0)
	cerr << "Nelze vypocitat NDVI" << endl;

    GDALClose(poDsTm3);
    GDALClose(poDsTm4);
    GDALClose(poDsNdvi);
    
    return 0;
}

Související články

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

Externí odkazy

• GDAL Data Model
• API Reference Documentation
• GDAL API Tutorial
• GDAL Driver Implementation Tutorial