-
Notifications
You must be signed in to change notification settings - Fork 0
/
calc_Indices_Ikonos.py
196 lines (164 loc) · 6.93 KB
/
calc_Indices_Ikonos.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Description:
This script uses a band set calculate the NDVI, the NDWI & the BuiltUp Index.
Afterwards the scripts eyport each single band plus the calculated indices.
Last the script merges the bands in a single band set (raster stack).
Input: wholeFilename.tif
Output: Sindle bands & indices as tif files &
Virtual & usual raster band set
Structure:
0. Functions
a) Export single bands exportSingleBand
b) NDVI calculator calcNDVI
c) BuiltUp Index calculator calcBuiltUp
d) NDWI Index Calculator calcNDWI
I. Import filename
II. Label postfix of the indices
III. Check if file exist
IV. Calculate NDVI
V. Calculate NDWI
VI. Calculate Built Up Index
VII. Label filenames for single bands & indices
VIII. Get one band and write it
IX. Label virtual & final raster
X. Build virtual raster
XI. Translate virtual raster to raster
Author: Mario Blersch
Date: February 2020
"""
# a) Function to Export single bands
def exportSingleBand(image_file, band_number, output_file):
# Create the file
with rasterio.open(image_file) as src:
band = src.read(band_number)
# Set spatial characteristics of the output object to mirror the input
kwargs = src.meta
kwargs.update(
dtype=rasterio.float32,
count = 1)
with rasterio.open(output_file, 'w', **kwargs) as dst:
dst.write_band(1, band.astype(rasterio.float32))
return 'exportSingleBand done'
# b) NDVI calculator
def calcNDVI(image_file, output_file):
# Formula source: https://wiki.orfeo-toolbox.org/index.php/Radiometric_Indices
# Load red and NIR bands - note all PlanetScope 4-band images have band order BGRN
with rasterio.open(image_file) as src:
band_red = src.read(3)
with rasterio.open(image_file) as src:
band_nir = src.read(4)
# Allow division by zero
numpy.seterr(divide='ignore', invalid='ignore')
# Calculate NDVI
ndvi = (band_nir.astype(float) - band_red.astype(float)) / (band_nir + band_red)
# Set spatial characteristics of the output object to mirror the input
kwargs = src.meta
kwargs.update(
dtype=rasterio.float32,
count = 1)
# Create the file
with rasterio.open(output_file, 'w', **kwargs) as dst:
dst.write_band(1, ndvi.astype(rasterio.float32))
return ndvi
# c) BuiltUp Index calculator
def calcBuiltUp(image_file, output_file):
# Source: https://wiki.orfeo-toolbox.org/index.php/ISU
# Load red and NIR bands - note all PlanetScope 4-band images have band order BGRN
with rasterio.open(image_file) as src:
band_red = src.read(3)
with rasterio.open(image_file) as src:
band_nir = src.read(4)
# Allow division by zero
numpy.seterr(divide='ignore', invalid='ignore')
# Calculate NDVI
builtup = 100 - 25 * band_nir.astype(float) / band_red.astype(float)
# Set spatial characteristics of the output object to mirror the input
kwargs = src.meta
kwargs.update(
dtype=rasterio.float32,
count = 1)
# Create the file
with rasterio.open(output_file, 'w', **kwargs) as dst:
dst.write_band(1, builtup.astype(rasterio.float32))
return builtup
# d) NDWI Index Calculator
def calcNDWI(image_file, output_file):
# Source: #https://wiki.orfeo-toolbox.org/index.php/NDWI_(Mc_Feeters,_1996)
# Load red and NIR bands - note all PlanetScope 4-band images have band order BGRN
with rasterio.open(image_file) as src:
band_green = src.read(2)
with rasterio.open(image_file) as src:
band_nir = src.read(4)
# Allow division by zero
numpy.seterr(divide='ignore', invalid='ignore')
# Calculate NDVI
ndwi = (band_green.astype(float) - band_nir.astype(float)) / (band_green.astype(float) + band_nir.astype(float))
# Set spatial characteristics of the output object to mirror the input
kwargs = src.meta
kwargs.update(
dtype=rasterio.float32,
count = 1)
# Create the file
with rasterio.open(output_file, 'w', **kwargs) as dst:
dst.write_band(1, ndwi.astype(rasterio.float32))
return ndwi
if __name__ == "__main__":
import rasterio
import numpy
import os.path
import os
# I. Import filename
img = '/home/geowazm/Documents/05_WS-201920/Master_Thesis/DATA/VHR_images/Ikonos/121-130/po_2661437_0000000/po_2661437_bgrn_0000000_2000_classified/po_2661437_bgrn_0000000_2000_img_clip.tif'
# II. Label postfix of the indices
label_ndvi = '_b5_ndvi.tif'
label_ndwi = '_b6_ndwi.tif'
label_builtup = '_b7_builtup.tif'
# III. Check if file exist
if os.path.isfile(img):
print ('{} -> exist'.format(img))
else:
print ("File not exist")
# IV. Calculate NDVI
ndvi = img[:-4] + '{}'.format(label_ndvi)
ndvi_result = calcNDVI(img, ndvi)
if ndvi_result.all is not None:
print('NDVI finished')
# V. Calculate NDWI
ndwi = img[:-4] + '{}'.format(label_ndwi)
ndwi_result = calcNDWI(img, ndwi)
if ndwi_result.all is not None:
print('NDWI finished')
# VI. Calculate Built Up Index
builtup = img[:-4] + '{}'.format(label_builtup)
builtup_result = calcBuiltUp(img, builtup)
if builtup_result.all is not None:
print('BuiltUp Index finished')
# VII. Label filenames for single bands & indices
blue = img[:-4] + '_b1_blue.tif'
green = img[:-4] + '_b2_green.tif'
red = img[:-4] + '_b3_red.tif'
nir = img[:-4] + '_b4_nir1.tif'
lbl_ndvi = img[:-4] + '_b5_ndvi.tif'
lbl_ndwi = img[:-4] + '_b6_ndwi.tif'
lbl_builtup = img[:-4] + '_b7_builtup.tif'
# VIII. Get one band and write it
#exportSingleBand(img, 1, cstl)
exportSingleBand(img, 1, blue)
exportSingleBand(img, 2, green)
exportSingleBand(img, 3, red)
exportSingleBand(img, 4, nir)
print('Single bands written')
# IX. Label virtual & final raster
label_vrt = '_vrt_stack.tif'
label_rst = '_rst_stack.tif'
vrt_raster = img[:-4] + '{}'.format(label_vrt)
raster_stack = img[:-4] + '{}'.format(label_rst)
# X. Build virtual raster
cmd_vrt = 'gdalbuildvrt -separate {} {} {} {} {} {} {} {}'.format(vrt_raster, blue, green, red, nir, lbl_ndvi, lbl_ndwi, lbl_builtup)
os.system(cmd_vrt)
# XI. Translate virtual raster to raster
cmd_final = 'gdal_translate {} {}'.format(vrt_raster, raster_stack)
os.system(cmd_final)
print('Bands merged')