GDAL / Programování / Rastrová data
Z FreeGIS portál
< 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.
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())
