-
Notifications
You must be signed in to change notification settings - Fork 0
/
sambuca_input_rrs.py
147 lines (102 loc) · 6.02 KB
/
sambuca_input_rrs.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
# -*- coding: utf-8 -*-
"""
Created on Mon Jan 30 16:40:19 2017
@author: Marco
"""
from os.path import join
import os.path
import rasterio
import sambuca as sb
import sambuca_core as sbc
#import nibabel as nib
#import snappy
import numpy as np
def sam_obs(base_path, Rrs = False):
if __name__=='sambuca_input_rrs':
#print (os.path.isfile('bioopti_data\\..\\sambuca\\reference\\wl_alos_data\\inputs\\WL_ALOS_R_0_sub120.img'))
#base_path = 'C:\\Users\\PCUSER\\sambuca_project\\input_data\\'
observed_rrs_base_path = base_path + 'image\\'
observed_rrs_raster_path = join(observed_rrs_base_path, 'WL_ALOS_R_0_sub120.img')
observed_rrs_filename='WL_ALOS_R_0_sub120.img'
sensor_filter_path = join(base_path, 'sensor_filters')
sensor_filter_name = 'ALOS'
nedr_path = join(base_path + 'nedr\\', 'WL_ALOS_NEDR_0_4bands.hdr')
observed_rrs_width = 0
observed_rrs_height = 0
observed_rrs = None
with rasterio.Env():
with rasterio.open(observed_rrs_raster_path) as src:
print('Observed rrs file: ', observed_rrs_raster_path)
print('Width, height: ', src.width, src.height)
print('crs: ', src.crs)
print('affine: ', src.affine)
print('num bands: ', src.count)
print('band indicies: ', src.indexes)
observed_rrs_width = src.width
observed_rrs_height = src.height
observed_rrs = src.read()
crs=src.crs
affine=src.affine
count=src.count
indexes=src.indexes
sensor_filters = sbc.load_sensor_filters(sensor_filter_path)
# We don't need to do this, but it lets us see the name of all loaded filters
sensor_filters.keys()
# retrieve the specified filter
sensor_filter = sensor_filters[sensor_filter_name]
# If Above surface remote sensing reflectance (Rrs) tag is true, convert to
# below surface (rrs) after Lee et al. (1999)
if Rrs == True:
observed_rrs = (2*observed_rrs)/((3*observed_rrs)+1)
nedr = sbc.load_spectral_library(nedr_path, validate=False)['wl_alos_nedr_0_4bands:33']
nedr
image_info={'observed_rrs_width':observed_rrs_width, 'observed_rrs_height':observed_rrs_height, 'crs':crs,
'affine':affine, 'count':count, 'indexes':indexes, 'nedr': nedr, 'sensor_filter':sensor_filter,
'base_path': base_path, 'observed_rrs_filename':observed_rrs_filename}
return observed_rrs, image_info
#return observed_rrs, observed_rrs_width, observed_rrs_height, nedr, sensor_filter, xstart, xend, ystart, yend, num_pixels, base_path, observed_rrs_filename
def sam_obs_s2():
if __name__=='sambuca_input_rrs':
print (os.path.isfile('bioopti_data\\..\\sambuca\\reference\\wl_alos_data\\inputs\\WL_ALOS_R_0_sub120.img'))
base_path = 'C:\\Users\\PCUSER\\sambuca_project\\bioopti_data\\'
observed_rrs_base_path = base_path + '..\\sambuca\\reference\\wl_alos_data\\inputs\\'
observed_rrs_raster_path = join(observed_rrs_base_path, 'S2_lampi_05021_atm_rrs')
observed_rrs_filename='S2_lampi_05021_atm_rrs'
sensor_filter_path = join(base_path, 'sensor_filters')
sensor_filter_name = 'ALOS'
nedr_path = join(observed_rrs_base_path, 'WL_ALOS_NEDR_0_4bands.hdr')
sensor_filter_path = join(base_path, 'sensor_filters')
sensor_filter_name = 'ALOS'
observed_rrs_width = 0
observed_rrs_height = 0
observed_rrs = None
# with rasterio.drivers():
with rasterio.open(observed_rrs_raster_path) as src:
print('Observed rrs file: ', observed_rrs_raster_path)
print('Width, height: ', src.width, src.height)
print('crs: ', src.crs)
print('affine: ', src.affine)
print('num bands: ', src.count)
print('band indicies: ', src.indexes)
observed_rrs_width = src.width
observed_rrs_height = src.height
observed_rrs_o = src.read()
crs=src.crs
affine=src.affine
count=src.count
indexes=src.indexes
observed_rrs=np.array([np.transpose(observed_rrs_o[0]), np.transpose(observed_rrs_o[2]), np.transpose(observed_rrs_o[3]), np.transpose(observed_rrs_o[6])])
sensor_filters = sbc.load_sensor_filters(sensor_filter_path)
# We don't need to do this, but it lets us see the name of all loaded filters
sensor_filters.keys()
# retrieve the specified filter
sensor_filter = sensor_filters[sensor_filter_name]
#Plot the sensor filter:
#plot_items.clear() #Python 3.3 and later only
nedr = sbc.load_spectral_library(nedr_path, validate=False)['wl_alos_nedr_0_4bands:33']
nedr
image_info={'observed_rrs_width':observed_rrs_width, 'observed_rrs_height': observed_rrs_height, 'crs': crs,
'affine': affine, 'count': count, 'indexes':indexes, 'nedr': nedr, 'sensor_filter':sensor_filter,
'base_path':base_path, 'observed_rrs_filename':observed_rrs_filename}
#return observed_rrs, observed_rrs_width, observed_rrs_height, nedr, sensor_filter, xstart, xend, ystart, yend, num_pixels, base_path, observed_rrs_filename
return observed_rrs, image_info