GDAL / Programování / Rastrová data

Z FreeGIS portál
Přejít na: navigace, hledání

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.


ndvi = \frac{(nir - vis)}{(nir + vis)}

kde nir je kanál družicové scény Landsat pro blízké infračervené spektrum, vis je viditelné červené spektrum, viz rozpis kanálů Landsat.

Testovací data ke stažení zde.

C++

Soubor ke stažení zde.

  1 #include <iostream>
  2 #include "gdal_priv.h"
  3 
  4 using std::cout;
  5 using std::cerr;
  6 using std::endl;
  7  
  8 GDALDataset *open_file(const char *filename)
  9 {
 10     GDALDataset *poDs;
 11     
 12     poDs = (GDALDataset *) GDALOpen(filename, GA_ReadOnly);
 13     if (poDs == NULL) {
 14         cerr << "Otevreni '" << filename << "' selhalo." << endl;
 15         return NULL;
 16     }
 17     
 18     return poDs;
 19 }
 20 
 21 int calculate_ndvi(GDALDataset *poDsNdvi, GDALDataset *poDsTm3, GDALDataset *poDsTm4)
 22 {
 23     GDALRasterBand *poBandNdvi, *poBandTm3, *poBandTm4;
 24     
 25     poBandNdvi = poDsNdvi->GetRasterBand(1);
 26     poBandTm3  = poDsTm3->GetRasterBand(1);
 27     poBandTm4  = poDsTm4->GetRasterBand(1);
 28 
 29     if (!poBandNdvi || !poBandTm3 || !poBandTm4)
 30 	return -1;
 31 
 32     unsigned char *pafScanlineTm3, *pafScanlineTm4;
 33     double *pafScanlineNdvi;
 34     int nXSize = poBandTm3->GetXSize();
 35     int nYSize = poBandTm3->GetYSize();
 36     if (nXSize != poBandTm4->GetXSize() ||
 37 	nYSize != poBandTm4->GetYSize())
 38 	return -1;
 39 
 40     pafScanlineTm3 = (unsigned char *) CPLMalloc(sizeof(unsigned char) * nXSize);
 41     pafScanlineTm4 = (unsigned char *) CPLMalloc(sizeof(unsigned char) * nXSize);
 42     pafScanlineNdvi = (double *) CPLMalloc(sizeof(double) * nXSize);
 43 
 44     for (int row = 0; row < nYSize; row++) {
 45 	poBandTm3->RasterIO(GF_Read, 0, row, nXSize, 1, 
 46 			    pafScanlineTm3, nXSize, 1, GDT_Byte, 
 47 			    0, 0);
 48 	poBandTm4->RasterIO(GF_Read, 0, row, nXSize, 1, 
 49 			    pafScanlineTm4, nXSize, 1, GDT_Byte, 
 50 			    0, 0);
 51 	
 52 	for (int col = 0; col < nXSize; col++) {
 53 	    pafScanlineNdvi[col] = double(pafScanlineTm4[col] - pafScanlineTm3[col]) /
 54 		(pafScanlineTm4[col] + pafScanlineTm3[col]);
 55 	}
 56 	poBandNdvi->RasterIO(GF_Write, 0, row, nXSize, 1, 
 57 			     pafScanlineNdvi, nXSize, 1, GDT_Float32, 
 58 			     0, 0);
 59     }
 60     
 61     CPLFree(pafScanlineTm3);
 62     CPLFree(pafScanlineTm4);
 63     CPLFree(pafScanlineNdvi);
 64     return 0;
 65 }
 66 
 67 int main(int argc, char **argv)
 68 {
 69     const char *filenameTm3, *filenameTm4;
 70      
 71     double trans[6];
 72     GDALDriver   *poDriver;
 73     GDALDataset  *poDsTm3, *poDsTm4, *poDsNdvi;
 74     
 75     if (argc != 3) {
 76 	cerr << "Pouziti: " << argv[0] << " tm3 tm4" << endl;
 77 	return 1;
 78     }
 79     
 80     filenameTm3 = argv[1];
 81     filenameTm4 = argv[2];
 82 
 83     // registrovat dostupne GDAL ovladace
 84     GDALAllRegister();
 85      
 86     // otevrit rastrove soubory 'tm3' a 'tm4' pro cteni
 87     poDsTm3 = open_file(filenameTm3);
 88     poDsTm4 = open_file(filenameTm4);
 89     if (poDsTm3 == NULL || poDsTm4 == NULL)
 90 	return 1;
 91 
 92     poDriver = poDsTm3->GetDriver();
 93     poDsNdvi = poDriver->Create("ndvi.tif", poDsTm3->GetRasterXSize(), poDsTm3->GetRasterYSize(),
 94 				1, GDT_Float32, NULL);
 95     poDsTm3->GetGeoTransform(trans);
 96     poDsNdvi->SetGeoTransform(trans);
 97     poDsNdvi->SetProjection(poDsTm3->GetProjectionRef());
 98     
 99     if (poDsNdvi == NULL)
100 	cerr << "Nelze vytvorit vystupni soubor 'ndvi.tif'" << endl;
101 
102     if (calculate_ndvi(poDsNdvi, poDsTm3, poDsTm4) != 0)
103 	cerr << "Nelze vypocitat NDVI" << endl;
104 
105     GDALClose(poDsTm3);
106     GDALClose(poDsTm4);
107     GDALClose(poDsNdvi);
108     
109     return 0;
110 }

Poznámka pro QtCreator

INCLUDEPATH += /usr/include/gdal
LIBS += -lgdal

Python

#!/usr/bin/env python

import numpy as np
from osgeo import gdal, osr

import sys

def open_tif(filename):
    ds = gdal.Open(filename)
    if ds is None:
        sys.exit("Nelze otevrit soubor {}".format(filename))

    band = ds.GetRasterBand(1)
    if band is None:
        sys.exit("Nelze nacist rastrovy kanal")

    band.SetNoDataValue(0)
    data = band.ReadAsArray()
    
    return ds, data

def create_tif(filename, ids, array):
    cols = array.shape[1]
    rows = array.shape[0]
    
    driver = gdal.GetDriverByName('GTiff')
    ds = driver.Create(filename, cols, rows, 1, gdal.GDT_Float32)

    band = ds.GetRasterBand(1)
    band.SetNoDataValue(0)
    band.WriteArray(array)

    ds.SetGeoTransform(ids.GetGeoTransform())

    srs = osr.SpatialReference()
    srs.ImportFromWkt(ids.GetProjectionRef())
    ds.SetProjection(srs.ExportToWkt())

    band.FlushCache()

    return ds

def compute_ndvi(nir, vis):
    # 16-bit Landsat 8 -> floats
    t = np.float32
    r = vis.astype(t)
    n = nir.astype(t) 
    
    # deleni nulou
    np.seterr(invalid='ignore')

    ndvi = (n - r) / (n + r)
    
    # zpatky do 16-bit
    ### ndvi = (ndvi + 1) * (2**15 - 1)
    
    return ndvi.astype(np.float32)

def main():
    if len(sys.argv) != 3:
        sys.exit("{} nir.tif vis.tif".format(sys.argv[0]))

    ds_nir, data_nir = open_tif(sys.argv[1])
    ds_vis, data_vis = open_tif(sys.argv[2])

    data_ndvi = compute_ndvi(data_nir, data_vis)

    ds_ndvi = create_tif("ndvi.tif", ds_nir, data_ndvi)

    return 0

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

Související články

Externí odkazy