-
Notifications
You must be signed in to change notification settings - Fork 2
/
dm.grass.sh
executable file
·197 lines (166 loc) · 6.23 KB
/
dm.grass.sh
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
197
#!/usr/bin/env bash
set -o errexit
set -o nounset
set -o pipefail
set -x
declare date=$1
declare infolder=$2
declare outfolder=$3
red='\033[0;31m'
orange='\033[0;33m'
green='\033[0;32m'
nc='\033[0m' # No Color
log_info() { echo -e "${green}[$(date --iso-8601=seconds)] [INFO] ${@}${nc}"; }
log_warn() { echo -e "${orange}[$(date --iso-8601=seconds)] [WARN] ${@}${nc}"; }
log_err() { echo -e "${red}[$(date --iso-8601=seconds)] [ERR] ${@}${nc}" 1>&2; }
trap ctrl_c INT # trap ctrl-c and call ctrl_c()
ctrl_c() { log_err "CTRL-C. Cleaning up"; }
trap err_exit ERR
err_exit() { log_err "CLEANUP HERE"; }
debug() { if [[ ${debug:-} == 1 ]]; then
log_warn "debug:"
echo $@
fi; }
[[ -d "${infolder}" ]] || (
log_err "${infolder} not found"
exit 1
)
mkdir -p "${outfolder}/${date}"
# Zoom to mask region
r.in.gdal input=mask.tif output=MASK --quiet
g.region raster=MASK
g.region zoom=MASK
g.region res=1000 -a
g.region -s # save as default region
# load all the data
yyyymmdd=${date:0:4}${date:5:2}${date:8:2}
scenes=$(
cd "${infolder}"
ls | grep -E "${yyyymmdd}T??????"
)
scene=$(echo ${scenes} | tr ' ' '\n' | head -n3 | tail -n1) # DEBUG
for scene in ${scenes}; do
g.mapset -c ${scene} --q
g.region -d --q
files=$(ls ${infolder}/${scene}/*.tif || true)
if [[ -z ${files} ]]; then
log_err "No files: ${scene}"
continue
fi
log_info "Importing rasters: ${scene}"
parallel -j 1 "r.external source={} output={/.} --q" ::: ${files}
log_info "Masking clouds in SZA raster"
r.grow input=SCDA_v20 output=SCDA_grow radius=-5 new=-1 --q # increase clouds by 5 pixels
# remove small clusters of isolated pixels
r.clump -d input=SCDA_grow output=SCDA_clump --q
# frink "(1000 m)^2 -> hectares" 100 hectares per pixel, so value=10000 -> 10 pixels
# this sometimes fails. Force success (||true) and check for failure on next line.
r.reclass.area -c input=SCDA_clump output=SCDA_area value=10000 mode=greater --q || true
[[ "" == $(g.list type=raster pattern=SCDA_area) ]] && r.mapcalc "SCDA_area = null()" --q
# SZA_CM is SZA but Cloud Masked: Invalid where buffered clouds over ice w/ valid SZA
r.mapcalc "SCDA_final = if((isnull(SCDA_area) && (MASK@PERMANENT == 220)) || (isnull(SCDA_v20) && (MASK@PERMANENT != 220)), null(), 1)" --q
r.mapcalc "SZA_CM = if(not(isnull(SZA)) & SCDA_final, 1, null())" --q
done
# Bands for which invalid pixels are kept (r_TOA and BT)
# bands=$(g.list type=raster mapset=* -m | grep -v PERMANENT | cut -d"@" -f1 | sort | uniq)
bands='BT_S7
BT_S8
BT_S9
height
mask
NDSI
O3
OAA
OZA
r_TOA_01
r_TOA_02
r_TOA_04
r_TOA_06
r_TOA_07
r_TOA_11
r_TOA_12
r_TOA_17
r_TOA_18
r_TOA_21
r_TOA_S1
r_TOA_S5
r_TOA_S5_rc
SAA
SCDA_area
SCDA_clump
SCDA_final
SCDA_grow
SCDA_v20
SZA
SZA_CM
WV'
g.mapset -c ${date} --quiet # create a new mapset for final product
r.mask raster=MASK@PERMANENT --o --q # mask to Greenland ice+land
# find the array index with the minimum SZA
# Array for indexing, list for using in GRASS
sza_arr=($(g.list -m type=raster pattern=SZA mapset=*))
sza_list=$(g.list -m type=raster pattern=SZA mapset=* separator=comma)
r.series input=${sza_list} method=min_raster output=sza_lut --o --q
# echo ${SZA_list} | tr ',' '\n' | cut -d@ -f2 > ${outfolder}/${date}/SZA_LUT.txt
# find the indices used. It is possible one scene is never used
sza_lut_idxs=$(r.stats --q -n -l sza_lut)
n_imgs=$(echo $sza_lut_idxs | wc -w)
# link sza_lut_idxs with scene IDs
output=$(
python <<END
import pandas as pd
import numpy as np
scenes_str = "${sza_list}"
scenes_sep = scenes_str.split(',')
scenes_number = [s.split('@')[-1] for s in scenes_sep]
idxs_ids = pd.DataFrame({'lut_index': np.arange(len(scenes_number)),
'scene': scenes_number})
print(idxs_ids)
END
)
# save idxs_ids in txt file
echo "${output}" >>${outfolder}/${date}/scenes_lut.txt
# generate a raster of nulls that we can then patch into
log_info "Initializing mosaic scenes..."
parallel -j 1 "r.mapcalc \"{} = null()\" --o --q" ::: ${bands}
### REFERENCE LOOP VERSION
# Patch each BAND based on the minimum SZA_LUT
# for b in $(echo $bands); do
# # this band in all of the sub-mapsets (with a T (timestamp) in the mapset name)
# b_arr=($(g.list type=raster pattern=${b} mapset=* | grep "@.*T"))
# for i in $sza_lut_idxs; do
# echo "patching ${b} from ${b_arr[${i}]} [$i]"
# r.mapcalc "${b} = if(sza_lut == ${i}, ${b_arr[${i}]}, ${b})" --o --q
# done
# done
# PARALLEL?
log_info "Patching bands based on minmum SZA_LUT"
doit() {
local lut=$1
local idx=$2
local band=$3
local b_arr=($(g.list type=raster pattern=${band} mapset=* | grep "@.*T"))
r.mapcalc "${band} = if((${lut} == ${idx}), ${b_arr[${idx}]}, ${band})" --o --q
}
export -f doit
parallel -j 1 doit {1} {2} {3} ::: sza_lut ::: ${sza_lut_idxs} ::: ${bands}
# diagnostics
r.series input=${sza_list} method=count output=num_scenes_cloudfree --q
mapset_list=$(g.mapsets --q -l separator=newline | grep T | tr '\n' ',' | sed 's/,*$//g')
raster_list=$(g.list type=raster pattern=r_TOA_01 mapset=${mapset_list} separator=comma)
r.series input=${raster_list} method=count output=num_scenes --q
bandsFloat32="$(g.list type=raster pattern="r_TOA_*") SZA SZA_CM SAA OZA OAA WV O3 NDSI BT_S7 BT_S8 BT_S9 r_TOA_S5 r_TOA_S5_rc r_TOA_S1 height SCDA_v20 SCDA_grow SCDA_clump SCDA_area SCDA_final"
bandsInt16="sza_lut num_scenes num_scenes_cloudfree"
log_info "Writing mosaics to disk..."
tifopts='type=Float32 createopt=COMPRESS=DEFLATE,PREDICTOR=2,TILED=YES --q --o'
parallel -j 1 "r.colors map={} color=grey --q" ::: ${bandsFloat32} # grayscale
parallel -j 1 "r.null map={} setnull=inf --q" ::: ${bandsFloat32} # set inf to null
parallel "r.out.gdal -m -c input={} output=${outfolder}/${date}/{}.tif ${tifopts}" ::: ${bandsFloat32}
tifopts='type=Int16 createopt=COMPRESS=DEFLATE,PREDICTOR=2,TILED=YES --q --o'
parallel "r.out.gdal -m -c input={} output=${outfolder}/${date}/{}.tif ${tifopts}" ::: ${bandsInt16}
# Generate some extra rasters
tifopts='type=Float32 createopt=COMPRESS=DEFLATE,PREDICTOR=2,TILED=YES --q --o'
r.mapcalc "ndbi = ( r_TOA_01 - r_TOA_21 ) / ( r_TOA_01 + r_TOA_21 )"
r.mapcalc "bba_emp = (r_TOA_01 + r_TOA_06 + r_TOA_17 + r_TOA_21) / 4.0 * 1.003 + 0.058"
r.out.gdal -f -m -c input=ndbi output=${outfolder}/${date}/NDBI.tif ${tifopts}
r.out.gdal -f -m -c input=bba_emp output=${outfolder}/${date}/BBA_emp.tif ${tifopts}